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ABSTRACT 

We  investigate  the  propagation  in  air  of  plane  finite -amplitude 
noise  waves.   The  reported  research  is  "both  theoretical  and  experimental, 
and  the  results  obtained  are  valid  both  "before  and  after  shock  formation. 
In  the  theoretical  portion  of  the  work,  a  new  model  called  modified 
weak-shock  theory  is  discussed.  In  this  model  the  effects  of  ordinary 
weak-shock  theory  and  Kirchhoff  tube  wall  attenuation  and  dispersion 
are  combined.  The  inclusion  of  tube  wall  effects  is  necessary  because 
the  experimental  work  is  done  in  a  pipe.  The  modified  weak-shock  theory 
is  implemented  in  a  computer  algorithm. 

Our  theoretical  model  is  novel  in  that  it  is  applicable  both 
before  and  after  shock  formation,  it  contains  the  effects  of  nonlinear 
distortion  as  well  as  tube  wall  attenuation  and  dispersion,  and  it  is 
applicable  for  input  waveforms  of  arbitrary  shape.   Previous  theories 
do  not  include  all  of  these  features. 

We  use  a  plane  wave  tube  in  our  experiments .  The  tube  is  96  ft 

long,  has  a  circular  cross  section,  and  is  terminated  in  such  a  way  that 

the  reflected  wave  is  ho   dB  lower  in  sound  pressure  level  than  the 

incident  wave.  The  acoustic  source  is  a  75  W  horn  driver  that  is  mounted 

at  one  end  of  the  pipe.  The  passband  of  this  source  is  700  Hz  to  2200  Hz 

and  it  can  generate  a  sound  pressure  level  of  160  dB  at  the  source.  We 

use  a  microphone  and  an  oscilloscope  to  observe  the  finite -amplitude  wave 

as  a  function  of  distance. 
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In  the  developmental  testing  of  this  new  model  -we  investigate 
the  propagation  of  sawtooth  waves  in  a  pipe.  Asymmetric  waveforms, 
reported  by  previous  researchers,  are  observed  and  compared  to  predictions 
made  using  the  new  theory.  We  obtain  excellent  agreement  between 
theory  and  experiment  over  the  frequency  range  500  to  3500  Hz  and  the 
sound  pressure  level  range  1^2  to  160  dB  re  0.0002  ubar. 

We  then  investigate  high  intensity  (155  to  160  dB)  noise 
pulses  and  we  observe  that  they  undergo  significant  distortion  as  a 
function  of  distance.   In  particular,  the  observed  pulses  simplify 
because  the  shock  fronts  merge  and  minor  waveform  irregularities  are 
fed  into  the  shocks .   The  waveform  at  the  source  is  used  as  an  input 
to  the  computer  algorithm.  Waveform  predictions,  obtained  via  the 
algorithm,  are  compared  to  the  observed  waveforms  and  excellent  agree- 
ment is  found. 

Computed  and  observed  noise  spectra  as  a  function  of  distance 
are  then  obtained.  The  spectra  show  significant  shaping  due  to  nonlinear 
effects.   The  spectra  obtained  at  locations  remote  from  the  source  are 
broader  in  frequency  content  than  the  source  spectrum  and  in  fact  the 
remote  spectra  actually  cross  over  the  source  spectrum.   In  this  work 
we  observe  a  high  frequency  crossover  at  about  5500  Hz  and  a  low  frequency 
crossover  at  about  350  Hz.   Hence,  we  arrive  at  the  very  important 
conclusion  that  nonlinear  propagation  can  cause  a  net  gain  in  energy  at 
frequencies  above  the  high  frequency  crossover  and  below  the  low 
frequency  crossover. 

We  conclude  this  work  with  a  discussion  of  the  physical  effects 

that  give  rise  to  the  observed  nonlinear  waveform  changes  and  spectral 

vi 


shaping.  We  also  include  a  discussion  of  the  role  of  dispersion  in  the 
propagation  of  finite -amplitude  waves  in  pipes. 
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CHAPTER  1    INTRODUCTION 

The  general  subject  of  this  dissertation  is  the  propagation 
of  plane,  finite -amplitude  sound  "waves  of  arbitrary  waveform.  We  are 
especially  interested  in  finite -amplitude  noise.   Our  approach  to  the 
subject  matter  is  both  experimental  and  theoretical,  and  the  results 
obtained  are  valid  both  before  and  after  shock  formation. 

We  restrict  our  attention  to  plane  waves,  even  though  most 
practical  problems  involve  nonplanar  waves.  The  rationale  behind  the 
plane  wave  restriction  is  quite  simple.  We  wish  to  check  experimentally 
our  theoretical  results  and  plane  waves  are  easily  generated  in  a  tube. 
A  secondary  benefit  is  a  simplification  of  the  mathematical  analysis. 
The  problem  solution  technique  that  we  use  is  implemented  in  a  computer 
algorithm  that  utilizes  existing  theoretical  results . 

Examples  of  high  intensity  sound  occur  in  many  industrial  and 
military  situations.   Ribner  (1969)  states  that  a  Saturn  V  rocket 
radiates  0.2  billion  acoustic  watts  and  generates  a  sound  pressure  level 
of  110  dB*  ten  miles  from  the  launch  site.  Webster  (1971)  reports 
sound  pressure  levels  as  high  as  1^0  to  150  dB  measured  on  the  flight 
deck  of  an  aircraft  carrier  50  ft  from  a  jet  airplane.   In  both  of  these 
examples  the  noise  source  is  broadband  in  frequency  content. 

The  preceding  examples  are  examples  of  high  intensity  sound, 
but  are  they  examples  of  the  type  of  sound  to  be  considered  in  this 
work?  To  answer  this  question,  we  must  carefully  define  what  we  mean 


Throughout  this  work,  sound  pressure  levels  are  referred  to 
0.0002  microbar.  -, 


when  we  use  the  term  finite -amplitude  distortion.   The  equations  that 
describe  the  propagation  of  sound  are  nonlinear  and,  in  general,  cannot 
be  solved.   If  the  magnitude  of  the  nonlinear  terms  in  the  equations 
is  small,  however,  the  equations  can  be  linearized  and  solved.   The 
applicability  of  linear  theory  is  amply  illustrated  by  the  fact  that  the 
use  of  linear  theory  gives  adequate  results  for  the  vast  majority  of 
acoustic  problems.  We  are,  however,  interested  in  acoustic  problems 
that  require  the  retention  of  the  nonlinear  terms .  We  further  restrict 
our  interest  to  those  problems  wherein  the  main  effect  of  nonlinearity  is 
a  cumulative  distortion  of  the  waveform.   To  be  sure,  there  are  important 
problems  in  which  the  primary  nonlinear  effects  are  local  rather  than 
cumulative;  but  in  the  wave  propagation  problems  we  shall  consider, 
noncumulative  nonlinear  effects  are  minor  in  comparison  to  the  cumulative 
ones.   Henceforth  in  this  dissertation  the  term  finite -amplitude  distor- 
tion should  be  considered  synonymous  with  the  term  cumulative  distortion. 
Now,  before  we  can  consider  the  question  posed  at  the  beginning  of  this 
paragraph,  we  must  identify  the  factors  that  affect  finite-amplitude 
distortion. 

The  magnitude  of  finite -amplitude  distortion  depends  on  three 
factors:   the  frequency  (or  frequencies)  of  the  disturbance,  the  local 
amplitude  of  the  disturbance,  and  the  propagation  distance  of  interest. 
The  sound  wave  is  distorted  over  each  wavelength  propagated  and  the 
total  amount  of  distortion  is  the  sum  of  the  distortion  that  occurs  in 
each  wavelength  (whence  the  term  "cumulative"  used  above).  Without 
going  into  detail,  inappropriate  in  an  introductory  chapter,  we  cannot 
say  with  certainty  whether  the  sound  waves  in  the  examples  discussed 


earlier  will  exhibit  appreciable  finite -amplitude  distortion.  We  can, 
however,  state  that  the  two  examples  appear  to  be  very  strong  candidates 
for  examples  of  the  type  of  sound  wave  in  which  we  are  interested. 

The  specific  objective  of  this  research  then  is  to  investigate 
the  propagation  of  plane,  finite -amplitude  noise.   The  noise  waveforms 
of  interest  are  broadband  in  frequency  content  and  intense  enough  that  the 
first  shocks  in  the  wave  develop  within  a  short  distance  from  the  noise 
source.   In  support  of  this  objective  we  develop  a  computer  algorithm 
that  can  be  used  to  predict  the  distortion  of  a  wave  whose  waveform  at 
the  source  is  arbitrary.   In  this  algorithm  we  calculate,  via  weak-shock 
theory,  the  distortion  suffered  by  a  wave  in  nonlinearly  propagating 
a  short  distance.  The  distorted  wave  is  then  corrected  for  tube  wall 
attenuation  and  dispersion  and  another  short  distortion  step  is  taken. 
The  waveform  for  any  value  of  propagation  distance  may  be  computed  in 
this  fashion.  Predictions  obtained  by  use  of  the  algorithm  are  then 
compared  to  experimental  measurements  made  using  a  plane  wave  tube. 


CHAPTER  2  REVIEW  OF  THE  LITERATURE 

1.  Introduction 

In  this  chapter  we  have  two  specific  objectives.   The  first 
is  to  review  the  development  of  the  theoretical  ideas  that  provide  the 
basic  underpinning  for  the  computer  methods  we  shall  develop.   The  second 
objective  is  to  briefly  consider  the  results  of  other  researchers  who 
have  worked  on  problems  similar  to  those  considered  here  and,  where 
applicable,  to  comment  on  how  their  ideas  have  helped  to  shape  and 
influence  our  own. 

2.  Background  of  the  Basic  Theory 

We  call  the  analytic  method  that  is  developed  in  this 
dissertation  modified  weak-shock  theory.   In  this  theory  three  factors 
are  combined  to  form  a  formal  method  for  considering  the  propagation  of 
a  finite -amplitude  sound  wave  in  a  tube.   These  factors  are  as  follows: 

a.  The  distortion  of  a  continuous  sound  wave  due  to  nonlinear 
effects.  The  formal  method  that  we  shall  use  to  calculate 
this  distortion  is  called  the  Earnshaw  (i860)  solution. 

b.  Once  the  sound  wave  is  no  longer  continuous  (i.e.,  shocks 
form)  the  Earnshaw  solution  must  be  supplemented  by  a 
relation  that  describes  the  shock  motion.   The  combination 
of  the  Earnshaw  solution  with  the  shock  motion  relation  is 
called  weak-shock  theory. 

c.  Since  our  experimental  work  will  be  done  using  a  plane  wave 
tube,  it  is  necessary  to  account  for  the  effects  of 
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attenuation  and  dispersion  due  to  the  tube  wall  boundary 
layer.  We  shall  modify  weak-shock  theory  in  an  iterative 
fashion  to  account  for  the  tube  wall  effects. 

Each  of  these  factors  will  now  be  briefly  described. 
a)   The  Earnshaw  Solution 

Poisson  (1808)  derived  the  classic  result  for  plane  sound 

waves  of  finite  amplitude, 


ij/     =   f[x+(a±\|>    )t]  ,  (2-1) 

X  X 


where 

\|/  =  the  velocity  potential, 

ty  =  the  derivative  of  i|/  with  respect  to  distance  x, 

X 

t  =  time,  and 

a  =  the  speed  of  sound  in  a  process  that  obeys  Boyle's  law. 
In  Eq.  (2-1)  the  upper  signs  apply  to  outgoing  waves,  the  lower  signs  to 
incoming  waves.  As  Blackstock  (1969)  points  out,  Poisson  missed  the 
importance  of  his  result.   Poisson  recognized  that  the  propagation 
velocity  dx/dt,  which  is  obtained  from  Eq.  (2-l),  depends  upon  \|/  .   But 
because  such  a  result  had  never  been  observed  experimentally,  he  concluded, 
"...  that  the  original  disturbance  is  transmitted  uniformly  and  with  a 
speed  equal  to  a  .,,."  We  now  recognize  that  most  acoustic  wave 
propagation  is  isentropic  and  that  a  is  properly  replaced  by  c,  where 


s 


and 

P  =  absolute  pressure, 

p  =  density,  and 

s  =  entropy. 
Thus  the  more  general  form  of  "Poisson's  solution"  is  Eq.  (2-l)  -with  a 
replaced  by  c. 

We  shall  be  concerned  here  -with  simple,  or  unidirectional, 
waves .  For  a  simple  wave  in  air  the  relation  between  the  local  sound 
speed  c  and  the  particle  velocity  u  is  [see,  for  example,  Courant  and 
Friedrichs  (19^8)] 

c  =  CQ  ±  1/2  (r-l)u    ,  (2-3) 

where 

c  =  the  small-signal  sound  speed,  and 

7  =  ratio  of  specific  heats. 
In  Eq.  (2-3)  the  upper  signs  are  used  when  considering  outgoing  waves 
and  the  lower  signs  when  considering  incoming  waves.  Substituting 
Eq.  (2-3)  into  Eq.  (2-1)  and  recognizing  that  \|r  =u,  we  obtain  the 
propagation  velocity 


§         =PutCo    ,  (2-M 

/  u= const 


where  fi=(y+l)/2   for  gases.   Equations  (2-3)  and  (2-^)  may  be  generalized 
to  apply  to  liquids  by  using  the  relation 

P  =  1  +  B/2A    ,  (2-5) 
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where  B  and  A  are  constants  that  depend  upon  the  liquid  being  considered 
[see,  for  example,  Beyer  (i960)].   The  dependence  of  the  propagation 
velocity  upon  the  particle  velocity,  as  described  "by  use  of  Eq.  (2-1+) , 
is  the  "basic  physical  phenomenon  causing  the  cumulative  distortion 
discussed  in  Chapter  1.  A  much  more  complete  discussion  of  Eq.  (2-4) 
appears  in  Chapter  3- 

Earnshaw  (i860)  considered  the  problem  of  the  lossless 
propagation  of  a  plane  wave  in  a  tube.   He  was  aware  of  the  variable 
propagation  velocity  implicit  in  Poisson's  result.  Earnshaw's  solution, 
which  we  postpone  giving  until  Chapter  3>  relates  the  shape  of  any 
continuous  wave  at  some  field  point  (x,t)  to  the  shape  of  that  wave  at 
the  source.   Hence,  the  Earnshaw  result  can  be  used  to  describe  the 
progressive  distortion  of  a  simple  plane  sound  wave. 
b)   Weak -Shock  Theory 

The  Earnshaw  solution  provides  an  adequate  method  for 
calculating  distortion  so  long  as  the  sound  wave  remains  continuous . 
In  general,  however,  the  Earnshaw  solution  leads  to  multivalued  wave- 
forms. Mult i value dness  signals  the  presence  of  shocks  in  the  waveform. 

An  analytic  method  called  weak-shock  theory,  which  is  basically 
an  extension  of  the  Earnshaw  solution,  enables  one  to  proceed  with  the 
calculation  of  waveform  distortion  when  shocks  are  present,  provided  the 
shocks  are  not  strong.   By  use  of  the  Rankine-Hugoniot  shock  relations, 
we  obtain  an  equation  that  describes  the  velocity  of  the  shock  front  in 
a  simple  wave.   The  propagation  velocity  and  attendant  distortion  in  the 
remainder  of  the  wave  is  calculated  using  the  Earnshaw  relation. 
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Several  researchers  contributed  to  the  development  of 
weak-shock  theory;  foremost  among  them  were  Landau  (19^5)>  Friedrichs 
(19^8),  and  Whitham  (1952  and  1956).   Blackstock  (1966)  stated  the 
theory  in  the  form  used  here.   The  adjective  "-weak"  refers  to  simplifying 
assumptions  concerning  the  shock  amplitude  made  in  the  theoretical 
development. 

c)   Tube  Wall  Effects 

The  experimental  observations  reported  here  were  made  using 
a  plane  wave  tube.  As  we  shall  see  later,  attenuation  and  dispersion 
caused  by  the  presence  of  the  tube  wall  boundary  layer  are  significant 
effects  and  must  be  considered  in  our  calculations.   The  only  attenuative 
effect  that  is  accounted  for  in  weak-shock  theory  is  dissipation  at  the 
shock.  We  shall  calculate  attenuation  and  dispersion  in  the  remaining 
portions  of  the  wave  using  the  classical  theory  first  derived  by 
Kirchhoff  (l868).   The  Kirchhoff  theory  is  a  "small-signal"  theory. 
Nonetheless,  it  can  be  combined  with  weak-shock  theory  to  form  what  we 
shall  call  modified  weak-shock  theory. 

3.   The  Contributions  of  Other  Researchers 

In  this  section  the  results  of  other  researchers  who  have 
worked  on  problems  similar  to  the  subject  of  this  study  are  briefly 
reviewed.   The  emphasis  in  this  review  is  on  the  relation  that  these 
other  studies  have  to  our  own.   The  section  is  organized  by  subjects 
as  follows : 

a.  General  studies  of  finite-amplitude  propagation  in  rigid  tubes. 

b.  Attenuation  of  sawtooth  waves  in  a  tube. 

c.  Computational  schemes. 


d.  Propagation  of  finite -amplitude  noise. 

e.  The  "Burgers  Turbulence"  model. 

a)   General  Studies  of  Finite -Amplitude  Propagation  in  Tubes 

The  problem  of  propagation  in  tubes  has  attracted  the  attention 
of  a  number  of  previous  researchers.  We  shall  discuss  in  chronological 
order  those  research  efforts  -which  we  consider  to  be  the  most  significant. 

Nearly  UO  years  ago  Thuras,  Jenkins,  and  O'Neil  (1935) 
investigated  the  propagation  of  an  initially  sinusoidal  tone  in  a  pipe. 
The  amplitude  and  frequency  of  the  tone  were  such  that  only  the  planar 
mode  was  excited  and  shocks  did  not  form  in  the  pipe.   Using  a  lossless 
wave  equation  in  Lagrangian  coordinates,  Thuras  et  al.  found  ordinary 
differential  equations  for  the  change  in  amplitude  of  the  second  and  third 
harmonics  of  the  fundamental  as  a  function  of  distance  from  the  source. 
The  differential  equation  for  the  second  harmonic  was  then  corrected 
for  the  attenuative  effect  (but  not  the  dispersive  effect)  of  the  tube 
wall  boundary  layer  by  postulating  that  both  the  fundamental  and  second 
harmonic  decrease  in  amplitude  at  a  rate  proportional  to  the  amount  of 
that  frequency  component  present.   The  differential  equation  proposed  was 


dPp      p  -205  x 

-^   =  KP  e   L     -  a   P     ,  (2-6) 

dx      o  2  2 


where 


Pp  =  the  second  harmonic  pressure  amplitude, 

P  =  the  fundamental  pressure  amplitude  at  x=0, 

05  and  o;  =  attenuation  coefficients  at  the  fundamental  and  second 

harmonic  frequencies  respectively,  and 
K  =  a  constant. 
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Predictions  obtained,  using  Eq.  (2-6)  were  compared  -with  experimentally 
observed  second  harmonic  levels .   For  the  propagation  distances 
considered,  the  second  harmonic  grew  linearly  with  distance.  Agreement 
between  theory  and  experiment  was  good  with  the  exception  of  a  constant 
3  dB  discrepancy.   Geertsen  (1951),  Blackstock  ( 1962b),  and  Cruikshank 
(1966)  reported  similar  experiments  in  which  good  agreement  between 
theory  and  experiment  was  found.   It  seems  likely  that  the  "missing 
3  dB"  of  Thuras,  Jenkins,  and  O'Neil  was  a  calibration  error. 

Using  a  perturbation  analysis,  Blackstock  ( 1962b)  demonstrated 
that  dispersion  as  well  as  attenuation  should  be  included  in  the 
analysis  of  the  second  harmonic.   Although  predictions  obtained  from  his 
expression  agreed  well  with  his  experimental  data,  the  effect  of  dispersion 
was  quite  small  for  the  conditions  of  his  experiment. 

Burns  (1966a)*  continued  Blackstock' s  perturbation  analysis, 
greatly  expanding  its  scope.   He  calculated  the  change  with  distance 
of  the  fundamental  through  the  fourth  harmonic.   His  theoretical  results 
include  the  effects  of  nonlinear  distortion  as  well  as  tube  wall  attenua- 
tion and  dispersion  but  are  limited  to  propagation  distances  that  preclude 
shock  formation.   Unfortunately,  however,  assumptions  made  to  simplify 
the  analysis  impose  a  further  limitation  on  the  propagation  distance 
over  which  his  theory  can  be  used.  A  secondary  result  of  these  assump- 
tions is  to  deemphasize  the  effect  of  dispersion  on  the  various  harmonic 
components.   However,  the  effects  of  dispersion  are  small  over  the 


For  a  shorter  and  more  accessible  reference,  see  Burns  (1966b). 
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limited  distances  that  were  considered.  Hence,  the  deemphasis  of  the 

role  of  dispersion  led  to  no  significant  errors.   If  much  longer 
propagation  distances  were  to  be  considered,  however,  Burns '  treatment 
of  dispersion  would  require  modification.  Burns  conducted  experiments 
meant  to  verify  his  theory.  The  propagation  distances  used  in  the 
experiments  were  limited  to  seven  wavelengths,  or  distances  less  than 
the  experimental  shock-formation  distance,  whichever  came  first. 
These  distances  were  far  too  short  for  dispersion  to  be  significant. 
The  agreement  between  theory  and  experiment  is  good  within  a  few  wave- 
lengths of  the  source,  but  only  fair  at  larger  distances. 

McKittrick,  Blackstock,  and  Wright  (1967)  presented  experimental 
results  that  clearly  showed  wave  propagation  phenomena  that  had  not  been 
previously  reported.  Oscillograms  of  shock  waves  formed  by  the  nonlinear 
distortion  of  sinusoids  showed  rounded  peaks  and  sharply  angled  troughs 
rather  than  the  expected  symmetric  sawtooth  shape.  The  authors  showed 
that  the  asymmetry  of  the  waveform  was  probably  due  to  the  influence  of 
dispersion,  caused  by  the  tube  wall  boundary  layer,  on  the  propagating 
wave. 

Pernet  and  Payne  (1969)*  also  considered  the  problem  of  the 
propagation  of  finite-amplitude  waves  in  pipes.  The  analytic  method 
employed  is  very  similar  to  that  used  by  Thuras,  Jenkins,  and  O'Weil. 


The  1969  report  appears  in  a  shorter  but  more  accessible  paper,  Pernet 
and  Payne  (l971a).  Unfortunately,  the  paper  does  not  contain  the 
discussion  of  dispersion  and  waveform  asymmetry  given  in  the  1969  report 
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Formulae  describing  the  change  in  amplitude  as  a  function  of  distance 
for  the  fundamental  through  the  seventh  harmonic  of  an  initially 
sinusoidal  -wave  were  derived.   Results  obtained  from  the  formulae  for 
the  fundamental  through  fourth  harmonic  were  experimentally  verified 
by  comparison  with  experimental  results  obtained  using  a  source  "wave  of 
moderate  (i.e.,  no  shock  formation)  intensity.  We  emphasize  here  that 
Pernet  and  Payne's  analytic  methods  are  only  applicable  for  distances 
prior  to  shock  formation. 

Pernet  and  Payne  mention  that  oscillograms  of  high  intensity 
waves  (e.g.,  160  dB)  show  the  type  of  asymmetry  described  earlier.   Like 
McKittrick  et  al.  they  conclude  that  dispersion  is  the  cause  of  the 
asymmetry.  Although  they  do  not  present  any  theory  applicable  to  waves 
of  high  intensity  they  do  present  some  experimental  measurements.   In 
Chapter  6  of  the  present  work  we  make  a  comparison  between  their  high 
intensity  experimental  data  and  predictions  based  on  modified  weak-shock 
theory.   Excellent  agreement  is  obtained.   Incidentally,  Burns'  theory 
does  not  agree  with  Pernet  and  Payne's  experimental  results,  nor  should 
it,  since  the  propagation  distances  considered  by  Pernet  and  Payne 
far  exceed  the  distance  limits  built  into  Burns'  theory. 

An  entirely  different  approach  from  those  presented  above  has 
been  proposed  by  Blackstock  (1972).   In  this  work  the  author  argues  that 
the  proper  wave  equation  describing  the  propagation  of  finite -amplitude 
waves  in  a  pipe  of  circular  cross  section  is 


du  fl  -2   on    1 
oi  "PCo   U  oV  =  R 


h 

nc 
o 
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(2-7) 
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where 

v  =  the  kinematic  viscosity, 

R  =  the  pipe  internal  radius, 

Pr  =  the  Prandtl  number,  and 

t'  =  t-x/c  ,  a  retarded  time. 
'  o 

The  right  side  of  Eq.  (2-7)  represents  the  effects  of  attenuation  and 
dispersion  caused  by  the  boundary  layer  and  the  u  du/dt'  term  represents 
the  effects  of  nonlinear  distortion.   It  is  not  clear  at  this  time 
whether  Eq.  (2-7)  could  be  used  alone  or  whether  use  of  the  Rankine- 
Hugoniot  relations  would  also  be  necessary.   In  either  case,  it  appears 
that  solving  Eq.  (2-7)  would  be  an  ideal  way  to  approach  the  problem 
of  finite -amplitude  propagation  in  a  pipe.  However,  the  analytic 
solution  of  Eq.  (2-7),  even  for  very  simple  source  functions,  has  not 
been  found.  An  alternative  approach  would  be  a  numerical  solution  of 
Eq.  (2-7)  obtained  by  means  of  a  computer.  This  path  has  not  been 
pursued  either. 

Coppens  (1971)  authored  the  last  paper  that  we  shall  consider 
in  this  section.   He  started  with  the  nonlinear  one -dimensional  wave 
equation  developed  by  Coppens  and  Sanders  (1968).  Coppens'  equation  is 
different  from  Eq.  (2-7)  in  form  but  similar  in  content.   The  equation 
is  solved  using  both  pertubation  techniques  and  Fourier  expansion.   The 
latter  method  is  amenable  to  computer  implementation  and  provides  the 
major  portion  of  Coppens'  analytic  results.   The  work  reported  is 
purely  theoretical  and  is  limited  to  the  region  before  shocks  form. 
Coppens'  computed  waveforms  clearly  exhibit  the  asymmetry  previously 
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reported  by  McKittrick  et  al.  (1967)  and  Pernet  and  Payne  (1969). 
Coppens  attributes  this  asymmetry  to  the  effects  of  dispersion. 

The  studies  Just  reviewed  are  of  value  here  for  the  theoretical 
and  experimental  background  they  provide,  rather  than  for  any  specific 
method  or  technique.  The  exception  to  this  statement  is  the  apparent 
necessity  of  including  the  effects  of  dispersion  in  our  analysis.  The 
necessity  might  be  regarded  as  questionable,  since  Thuras,  Jenkins, 
and  O'Neil  (1935),  Geersten  (1951),  and  Pernet  and  Payne  (1969)  ignored 
dispersion  and  nonetheless  got  good  agreement  between  theory  and 
experiment.  We  shall  include  dispersion  in  our  analysis,  and  in  our 
concluding  remarks  in  Chapter  7  indicate  why  dispersion  is  important 
in  some  analyses  but  may  be  ignored  in  others. 

b)   Attenuation  of  Sawtooth  Waves  in  a  Tube 

A  major  portion  of  the  study  of  finite-amplitude  sound  in 
pipes  has  been  directed  at  the  specific  problem  of  the  attenuation 
of  sawtooth-shaped  plane  waves .  The  various  theoretical  analyses  of 
this  problem  have  concluded  with  predictions  of  the  attenuation  per 
wavelength  propagated.   In  the  several  analyses  that  we  shall  review 
all  of  the  authors  took  account  of  attenuation  due  to  irreversible 
compression  that  the  fluid  medium  experiences  during  the  passage  of 
a  shock.  Additionally  all  of  the  authors  save  two  took  account  of 
attenuation  due  to  the  boundary  layer.  The  dependent  variable  commonly 
used  in  studies  of  sawtooth  wave  decay  is  [see  Rudnick  (1952)  or 
Medwin  (I96I+)] 

5  -A— '         >  (2-8) 

1 


15 


where 

P  =  the  total  pressure  at  the  wave  peak,  and 

P  =  the  total  pressure  at  the  wave  trough. 
The  value  of  6  at  the  space  origin  is  denoted  by  the  symbol  5  .   The 
independent  variable  is  the  dimensionless  distance  x/X,   where  X   is  the 
wavelength  of  the  sawtooth.   It  has  been  common  to  present  the  experi- 
mental data  by  plotting  5   versus  x/X.      Then  a  convenient  quantity  that 
describes  the  relationship  between  5  "  and  x/X   is  the  slope  S,  given  by 

Rudnick  (1952)  derived  the  value,  which  we  shall  call  S  , 

\-y-^r     ■  (2-10) 

Equation  (2-10)  may  also  be  derived  using  weak-shock  theory.   Fay  (1956) 
proposed  a  new  value  of  S,  viz., 


S2  =  ^4    ,  (2-U) 


but  later  withdrew  this  result,  and  affirmed  the  value  of  S  derived  by 
Rudnick,  S  [see  Fay  (1962)]. 

Westervelt  (i960)  concluded  that  the  proper  value  of  S  is 
given  by 

S,  =  ^±  .  (2-12) 

0  27 
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Temkin  (1969)  included  the  attenuative  effects  of  the  tube  wall  directly 
in  his  value  of  S  and  obtained 


s^^  +  HM)     >  <2-»: 


where  a=the  tube  wall  attenuation  coefficient  for  the  fundamental. 
Rudnick  (1952),  on  the  other  hand,  included  tube  wall  attenuation  in  his 
analysis  by  adjusting  his  experimental  data  in  such  a  fashion  that  tube 
wall  attenuation  losses  are  removed.   This  adjustment  is  accomplished 
using  a  correction  term  based  upon  the  attenuation  coefficient  of  the 
fundamental.   The  correction  term  is  commonly  referred  to  as  "Rudnick' s 
zero-order  correction."  Westervelt  ignored  tube  wall  attenuation 
entirely  in  his  analysis.   None  of  the  above  researchers  considered  the 
effects  of  dispersion. 

Both  Temkin  and  Rudnick  gave  the  limits  on  the  applicability 
of  their  analyses.  Temkin  stated  that  his  result  S,  is  valid  for  values 
of  6  <0.1,  whereas  Rudnick  said  simply  that  his  expression  S  is  valid 
for  "5  of  order  0.1."  Westervelt  did  not  mention  a  limit  to  the  applica- 
bility of  S.,.   The  maximum  value  of  8  is  an  important  limit  in  our 
analysis  because  it  can  be  used  to  identify  the  largest  sound  pressure 
level  that  can  be  accurately  investigated  using  weak-shock  theory.   Our 
analysis  of  this  limit  appears  in  Appendix  E. 

Hence,  at  present,  there  are  at  least  three  distinct  values 
of  S(i.e.,  S  ,  S  ,  and  S.  )  supported  by  theoretical  arguments. 
Unfortunately,  conclusive  experimental  confirmation  does  not  exist  for 
any  of  these  values . 
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Using  a  10  in.  plane  wave  tube,  Rudnick  (1952)  "with  experimental 

assistance  from  Leonard  and  Wilson  (1952)  investigated  the  attenuation 

of  intense  (O.O76  to  O.25  atm  peak  acoustic  pressure),  low  frequency 

(30  to  200  Hz)  sawtooth  waves.   The  peak  pressures  corresponded  to  a 

range  of  5  of 
0      o 

O.167  <  8  <  O.667    •  (2-1^) 

o 

Rudnick  found  that  the  observed  value  of  S  was  less  than  that  given 
by  Eq.  (2-10).   In  other  words,  the  amplitude  of  the  observed  sawtooth 
waves  decreased  at  a  slower  rate  than  predicted.  Adjustment  of  the 
data  using  the  zero-order  correction  was  accomplished  but  the  effect 
was  slight,  since  for  Rudnick 's  experimental  conditions  the  correction 
terms  are  very  small.   Rudnick  formed  ratios  of  the  observed  values  of 
S  and  S  .  With  a  single  exception,  these  ratios  fell  between  0.5  and 
1.0,  with  the  majority  near  0.7.  It  is  noteworthy  that  the  ratios 
appeared  to  tend  towards  1.0  as  the  value  of  6  decreased  toward  0.1, 
thereby  suggesting  that  the  value  S  is  valid  only  for  values  of  8  §0.1. 
Rudnick' s  experimental  methods  were  extended  to  other  frequency  and 
peak  pressure  ranges,  other  gases,  and  other  pipe  sizes  by  Werth  and 
Delsasso  ( 195*0  an<i  Medwin  ( 196*0 .   In  general  the  results  reported 
are  similar  to  Rudnick' s.  That  is,  the  observed  values  of  S  were  less 
than  S  even  after  correction  for  tube  wall  losses.  Although  Werth 
et  al.  and  Medwin  extended  the  range  of  measured  results,  with  a  single 
exception  all  of  their  reported  data  was  gathered  using  sawtooth  waves 
for  which  8  X).l.   The  exception  is  a  portion  of  the  data  reported  by 
Medwin  (196U)  in  his  Fig.  K.      Temkin  (1969)  compared  his  predicted 
value  S,  to  this  data  and  obtained  favorable  agreement.  Medwin  also 
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reported  that  some  results  obtained  in  air  and  argon  confirm  Westervelt's 

value  S  ,  but  measurements  obtained  using  other  gases  fail  to  confirm 
3 

either  Sn  or  S_. 
1     3 

A  piece  of  evidence  in  support  of  validity  of  the  value  S 
has  been  reported  by  Shin  (1963).   Shin  recast  the  Rudnick  theory  in 
spherical  coordinates,  calculated  the  rate  of  change  of  6   with  distance, 
and  adjusted  this  calculation  for  spherical  spreading.   The  resultant 
theoretical  value  S  was  then  compared  to  experimental  values  of  S 
observed  in  an  anechoic  medium.   Shin  obtained  excellent  agreement 
between  theory  and  experiment.   The  range  of  values  of  6  ,  measured 
at  the  point  of  sawtooth  wave  formation,  was 

0.027  <  5  <  O.OU7    .  (2-15) 

In  Appendix  E  we  have  attempted  to  improve  on  Rudnick' s 

computation  of  the  quantity  d(l/5)/d(x/\)  by  using  the  exact  expression 

for  the  entropy  change  across  any  shock. *  On  the  basis  of  our  improvement 

of  Rudnick' s  result  we  conclude  that  the  weak-shock  value  S  is  accurate 

only  if  5  <0.1.   For  cases  where  5  X).l,  terms  ignored  in  Rudnick ' s 
o  o 

original  derivation  become  significant  and  must  be  considered.   Implicit 

in  both  Rudnick' s  derivation  and  ours  is  the  assumption  of  simple  wave 

flow.   However,  for  values  of  5  X).l  reflections  from  the  shock  may 

o 

cause  one  to  challenge  seriously  the  validity  of  this  assumption. 
Therefore,  our  result  should  not  be  considered  a  complete  correction  of 


Rudnick  used  a  one-term  approximation,  namely  the  first  term  in 
Eq.  (3-17). 
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Rudnick's  theory  but  rather  a  working  limit  on  the  range  of  applicability 
of  weak-shock  theory.  We  have,  nonetheless,  compared  both  our  three- 
term  result  and  numerical  integration  result  (see  Appendix  E)  with 
Rudnick's  (1952)  experimental  data  and  obtained  generally  good  agreement. 
c)   Computational  Schemes 

The  theory  developed  in  this  dissertation  has  been  implemented 
by  means  of  a  computer  algorithm.   Hence,  it  is  appropriate  that  we 
review  briefly  earlier  computational  schemes.   Fox  and  Wallace  ( 195*0 
investigated  the  distortion  of  an  initially  sinusoidal  wave  in  a 
dissipative  (but  not  dispersive)  medium  for  propagation  distances 
prior  to  the  formation  of  shocks .   Their  analytic  method  was  an  iterative 
numerical  scheme  wherein  a  separate  difference  equation  was  used  for 
each  harmonic.   The  significance  of  their  work  in  relation  to  the  present 
work  is  that  we  have  also  used  an  iterative  procedure  wherein  we  correct 
the  computed  waveform  for  attenuation  and  dispersion  after  each  incre- 
mental propagation  step.  A  procedure  not  limited  to  distances  that 
preclude  shock  formation  was  developed  by  Cook  (1962).  A  parametric 
integration  method  similar  to  that  used  to  generate  the  Bessel-Fubini 
solution  [see  Blackstock  ( 1962a)]  was  used.  The  basic  procedure  was 
iterative  with  absorption  (but  not  dispersion)  introduced  after  each 
incremental  distance  step.  An  extension  of  Cook's  method  to  nonsinu- 
soidal  deterministic  waves  was  formulated  by  Van  Buren  and  Breazeale 
(1968)  in  their  study  of  the  reflection  of  finite -amplitude  waves.   In 
both  these  algorithms  dissipative  effects  such  as  viscosity  are  relied 
upon  to  prevent  the  waveform  from  becoming  multivalued.   There  is  no 
guarantee  that  pipe  wall  effects  are  sufficient  by  themselves  to  prevent 
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multivaluedness  during  propagation  in  a  pipe.   Hence,  Cook's  method 
"would  have  to  be  significantly  modified  to  be  applicable  to  our  problem. 
A  second  difficulty  -with  Cook's  method  is  that  it  does  not  contain  any 
provision  for  relative  motion  of  shocks  -within  the  -waveform.  We  -will 
see  in  Chapter  3  that  in  an  arbitrarily  shaped  wave  shocks  exhibit 
relative  motion  -with  respect  to  one  another  and  with  respect  to  the 
continuous  portion  of  the  -wave. 

A  computer  model  -wherein  relative  motion  of  the  shocks  is 
allowed  has  been  formulated  by  Hawkings  (l970)«  An  initial  waveform 
that  is  sawtoothed  with  one  or  more  of  the  sawtooth  peaks  larger  than 
the  rest  is  postulated.  The  model  is  based  on  plane  wave  weak-shock 
theory,  and  the  waveform  is  simplified  by  representing  it  in  straight 
line  segments.   Of  all  the  computer  models  of  finite-amplitude  propaga- 
tion of  which  we  are  aware,  Hawkings'  model  is  the  most  similar  to 
ours.   The  major  differences  between  our  model  and  Hawkings'  are  the 
straight  line  representation,  which  we  do  not  require,  and  our  inclusion 
of  explicit  attenuation  and  dispersion  effects,  which  are  ignored  by 
Hawkings . 

d)   The  Propagation  of  Finite -Amplitude  Noise 

There  have  been  only  a  few  reported  research  efforts  devoted 
to  the  investigation  of  finite-amplitude  noise  propagation.   Burnett 
and  Ackerman  (i960)  considered  l/2  octave  wide  bands  of  noise,  and  claim 
validity  for  their  theoretical  analysis  for  propagation  distances 
insufficient  for  the  formation  of  shocks .   The  authors  arbitrarily 
imposed  a  maximum  propagation  limit  of  3  meters.  Using  this  limit, 
Burnett  and  Ackerman  were  able  to  define  frequency-source  amplitude 
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regions  in  which  l)  small-signal  approximations  are  valid,  2)  difference 
frequency  signals  are  important,  3)  sum  and  difference  frequency  signals 
are  important,  k)   the  proposed  theory  failed.   Their  analysis  may  be 
useful  for  very  short  ranges  and  source  conditions  that  preclude  shock 
formation  within  these  ranges  but  is  not  useful  for  the  present  study. 
The  results  of  Burnett  and  Ackerman's  theoretical  study  were  investigated 
experimentally  by  Oda  and  Ackerman  (i960).   Oda  and  Ackerman  claim  to  have 
verified  Burnett  and  Ackerman's  theory  in  regions  1  and  2. 

In  a  portion  of  the  previously  mentioned  work  by  Pernet  and 
Payne  (1969),  and  in  a  later  paper  by  the  same  authors  [ Pernet  and 
Payne  (l971t>)],  finite -amplitude  noise  propagation  in  tubes  was  considered. 
In  both  papers  the  authors  obtained  theoretical  results  and  described 
verifying  experiments .   In  the  first  paper  the  authors  compared  pure 
tones  to  narrow  bands  of  noise  centered  at  the  frequency  of  the  tone. 
The  energy  content  of  the  two  signals  was  equal.   The  harmonics  of  the 
pure  tone  were  calculated.   Since  the  amplitude  peaks  in  a  narrowband 
noise  signal  have  a  Rayleigh  distribution  [see,  for  example,  Horton  (1969)], 
Pernet  and  Payne  were  able  to  modify  their  theory  and  calculate  the 
statistical  amplitude  distribution  of  the  second  through  seventh  harmonics 
of  the  finite-amplitude  primary  band.   Over  the  distances  considered, 
it  was  shown  that  noise  harmonic  bands  grew  more  quickly  with  propaga- 
tion distance  than  did  pure  tone  harmonics.   The  theoretical  predictions 
were  experimentally  verified.   As  in  the  analysis  of  pure  tones,  Pernet 
and  Payne's  noise  theory  is  limited  to  distances  that  preclude  shock 
formation.   In  the  second  paper  the  authors  utilized  the  fact  that  the 
amplitude  of  the  second  harmonic  is  proportional  to  the  squared  amplitude 
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of  the  fundamental.   This  quadratic  relation  was  analyzed  using  the 
mathematics  of  the  square-law  detector.   The  derived  expression  for  the 
second  harmonic  is  valid  for  distances  before  the  formation  of  shocks. 
The  theory  was  experimentally  verified.   Since  we  want  to  consider 
broadband  noise  and  are  interested  in  a  solution  valid  both  before  and 
after  shock  formation,  we  require  more  general  methods  than  those  proposed 
by  Pernet  and  Payne.   Therefore,  we  shall  not  make  specific  use  of  their 
analytic  results. 

Kuznetsov  (197O)  studied  the  spectra  of  high  intensity  noise. 
The  basis  for  his  theoretical  results  was  the  Fourier  transform  of  the 
Earnshaw  relation.   For  large  values  of  frequency  Kuznetsov  computed 
an  approximate  form  of  the  spectrum.   Although  Kuznetsov  does  not  mention 
it,  his  result  is  also  applicable  for  large  values  of  distance.   The 
results  are  somewhat  perplexing,  however,  because  Kuznetsov' s  theory  is 
only  applicable  for  distances  at  which  there  is  no  shock  formation. 

There  are,  of  course,  other  studies  wherein  the  authors  deal 
with  finite-amplitude  noise  [see,  for  example,  Jones  (1971)  or  Pickett 
(l97l)]  but  no  researchers  that  we  are  aware  of,  other  than  those 
commented  on  above,  have  dealt  directly  with  the  analysis  of  the  propaga- 
tion of  intense  random  noise. 

e)   The  "Burgers  Turbulence"  Model 

Burgers  (19^8)  proposed  a  mathematical  model  to  describe 
turbulence.  Using  this  model,  we  may  consider  the  development  of  some 
initial  value  turbulent  velocity  field  as  functions  of  a  space 


23 


coordinate  x  and  time  t.   The  ranges  of  the  independent  variables  in  the 
case  in  -which  we  are  interested  are  0  -»  °°.   The  equation  is 


Uj.  +  uu  =  vu      , 
t     x     xx 


(2-16) 


■where  v  represents  the  frictional  effects.   In  a  later  paper  by  Burgers 
(1972)  solutions  of  Eq.  (2-l6)  for  which  the  initial  velocity  field 
u(x, 0)  is  random  and  v  vanishingly  small  are  considered.  Under  these 
circumstances  the  solution  of  Eq.  (2-16)  for  large  values  of  t  assumes 
the  form  of  a  series  of  discontinuities  in  velocity  separated  by 
straight  line  sections  "whose  slope  is  t   .We  shall  refer  to  such  a 
solution  as  a  "random  sawtooth"  [see  Fig.  (2-l)]. 


Figure  2- 1 

A  sawtooth  pattern  as  the  asymptotic  solution 

of  u(x,  t)  for  large  t  and  v  vanishingly  small. 
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Although  Burgers'  results  were  originally  meant  to  apply  to 
turbulence,  Morfey  (1973)*  an(i  Pelinovskii  and  Fridman  (1973)  have 
noted  that  they  apply  equally  well  to  acoustic  noise  •waveforms  because 
Burgers'  equation  is  also  a  model  for  finite -amplitude  wave  propagation. 
In  the  papers  to  which  we  shall  refer  in  the  remainder  of  this  section 
the  authors  were  interested  in  turbulence  but  their  results,  because 
they  are  based  upon  Burgers'  equation,  may  also  be  construed  as  applying 
to  finite -amplitude  noise. 

Kraichnan  (1968)  noted  that  for  the  case  of  very  small  values 
of  v  (Kraichnan  used  the  term  "very  large  Reynolds  number")  that  the 
random  sawtooth  pattern  changes  with  time  because  discontinuities 
(shocks)  coalesce.   In  addition,  fluid  elements  travel  at  constant 
velocity  until  they  meet  a  shock  whereupon  they  disappear  into  the 
shock.   This  behavior  indicates  that  the  turbulent  field  should  simplify 
with  time.   Because  of  the  shock  coalescence,  the  spectrum  should 
show  a  shift  of  energy  to  lower  wave  numbers  [the  same  prediction  was 
made  earlier  by  Ffowcs  Williams  (1965)].  We  will  see  in  Chapter  6  that 
the  simplification  of  the  particle  velocity  field  is  as  described  by 
Kraichnan.   In  addition,  we  will  see  that  there  is  cascading  of  energy 
to  both  the  high  and  low  ends  of  the  noise  spectrum,  the  downward 
cascading  being  as  described  by  Kraichnan. 

The  random  sawtooth  waveforms  predicted  by  Burgers,  and  the 
merging  of  shocks,  waveform  simplification,  and  downward  shifting  of 
energy  predicted  by  Kraichnan  will  be  confirmed  later  in  this  work. 
However,  our  analytical  tool  is  weak-shock  theory,  not  Burgers'  equation. 
The  fact  that  use  of  weak-shock  theory  will  produce  results  equivalent 
to  asymptotic  solutions  of  Burgers'  equation  is  not  surprising,  since 
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it  is  known  that  the  Burgers '  equation  approach  and  -weak-shock  theory  are 
equivalent  for  the  case  of  v  -»  0  [see,  for  example,  Lighthill  (1956)]. 

Burgers' results  have  been  extended  by  Tatsumi  and  Kida  (1972) 
in  a  paper  devoted  primarily  to  the  statistical  properties  of  Burgers' 
model  of  turbulence.   Two  results  produced  by  Tatsumi  and  Kida  are 
equations  for  the  number  of  shocks  and  the  mean  interval  between  shocks 
as  functions  of  time.   These  relations  are 

N  .  VJ^j  (2-18) 

and 

<  =  <o(rJ     •  <2-19) 

where 

I     =   the  mean  interval  between  shocks, 

I     =   the  mean  interval  between  shocks  at  t=t  , 
o  o 

N  =  the  number  of  shocks, 

N  =  the  number  of  shocks  at  t=t  ,  and 
o  o 

i  =  a  constant,  £/t   is  the  number  of  shock  collisions  per 
unit  time. 
Via  statistical  arguments,  Tatsumi  and  Kida  show  that  1=1/ 2  corresponds 
to  the  condition  of  maximum  randomness  and  is  therefore  the  value  most 
likely  to  occur  in  real  turbulence.  A  numerical  calculation  of  the 
energy  of  Burgers  turbulence  was  carried  out  by  Crow  and  Canavan  (1970)* 
Tatsumi  and  Kida  extrapolated  the  value  of  I  from  the  Crow  and  Canavan 
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data  for  the  condition  of  very  large  Reynolds  number  and  obtained  the 
value  l/2. 

Burgers'  equation  is  convenient  for  initial  value  problems. 
For  boundary  value  problems  the  roles  of  "space"  and  "time"  in  the 
Burgers'  equation  are  reversed  [see,  for  example,  Blackstock  ( 1972)1. 
Recasting  Eq.  (2-l8)  so  that  it  applies  in  a  boundary  value  problem, 
we  obtain 


'x 

N  =  N  (  — )     .  (2-20) 

o\x 


where  x  is  a  reference  distance.  Utilizing  some  of  our  experimental 
results  obtained  using  noise  waves,  we  attempted  to  verify  Tatsumi  and 
Kida's  asymptotic  value  of  |.   The  results  of  this  attempt  are  reported 
in  Chapter  6. 

As  we  shall  see  in  Chapter  6,  we  are  unable  to  verify  Tatsumi 
and  Kida's  predicted  value  of  i  by  using  waves  of  initial  strength 
suitable  for  analysis  via  weak-shock  theory.   In  addition  the  asymptotic 
results  that  form  the  basis  of  the  work  of  Burgers,  Tatsumi  and  Kida, 
and  Kraichnan  appear  to  require  very  strong  initial  disturbances. 
Hence,  in  light  of  the  experimental  results  in  Chapter  6  and  the  analysis 
in  Appendix  E,  we  question  whether  weak-shock  theory  can  be  applied  to 
disturbances  of  the  strengths  required  to  obtain  the  asymptotic  results 
discussed  here. 


CHAPTER  3   THEORY 

1.   Introduction 

The  basic  -weak-shock  theory  used  in  this  dissertation  is 
■well  known  and  has  been  available  in  the  literature  for  many  years 
[see,  for  example,  Whitham  (1956),  or  Blackstock  (1966)  and  (1972)]. 
The  basic  equations  of  the  theory  are  briefly  derived  and  explained  in 
this  chapter  in  order  that  the  present  "work  be  complete  and  clear.   The 
application  of  -weak-shock  theory  is  then  demonstrated  in  two  example 
problems.  Results  similar  to  those  found  for  the  second  example  problem 
are  then  presented  for  a  generalized  plane  shock  propagation  problem. 
The  results  of  the  generalized  problem  illustrate  that  analytical 
solutions  of  even  fairly  simple  problems  are  quite  complex.  This 
complexity  suggests  that  a  computer  method  may  be  more  useful  than  an 
analytic  method  for  the  analysis  of  complicated  waveforms.  The  results 
of  the  generalized  problem  serve  a  second  purpose.  In  Chapter  2  we 
discussed  some  of  the  results  that  have  been  obtained  using  the  Burgers 
model  of  turbulence.   In  this  chapter  we  shall  show  that  equivalent 
results  may  be  obtained  from  the  asymptotic  forms  of  the  results  of 
the  generalized  problem. 

Our  theory  is  cast  in  terms  of  particle  velocity,  but  we 
actually  measure  pressure.  In  Section  h   of  this  chapter  we  develop 
the  necessary  relations  between  pressure  and  particle  velocity.  Next, 
we  consider  in  greater  detail  the  effects  of  the  tube  wall.  The  chapter 
is  concluded  with  a  discussion  of  modified  weak-shock  theory. 
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2.   Weak-Shock  Theory* 

a)   The  Earnshaw  Solution 

Consider  a  continuous  plane  sound  -wave  propagating  in  an 
otherwise  quiescent  medium.  We  assume  that  simple  wave  flow  exists, 
that  the  thermodynamic  processes  involved  are  isentropic,  and  that  the 
medium  is  a  perfect  gas . 

At  this  point  it  is  convenient  to  introduce  the  term  "-wavelet." 
A  wavelet  is  a  specific  point  on  a  waveform;  associated  with  this  point 
is  a  specific  particle  velocity  u.   If  more  than  one  point  on  the  wave 
has  the  same  particle  velocity,  each  point  is  still  a  different  wavelet. 

The  propagation  velocity  of  a  wavelet  of  particle  velocity  u 
is  given  by  Eq.  (2-4).   Restricting  Eq.  (2-k)   to  apply  to  only  outgoing 
waves,  we  obtain 

=  cq  +  pu    .  (3-1) 

u= const 

Equation  (3-l)  may  be  integrated  to  yield  the  Earnshaw  solution  once 

we  have  specified  appropriate  boundary  conditions.*-*  Let  us  assume  that 

a  plane  sound  wave  of  arbitrary  shape  is  being  observed  in  a  tube  by 

use  of  a  microphone  sensitive  to  particle  velocity.  We  will  consider  the 

microphone  location  the  origin  x-0  and  assume  that  the  wave  propagation 

is  lossless. 


Blackstock  (1972)  considers  waves  of  SPL  of  I'jk   dB  or  below  as  weak. 
Temkin  (1969),  in  a  more  detailed  analysis,  estimates  the  upper  limit 
of  weak-shock  theory  as  165  dB  or  below.  We  will  restrict  our  analysis 
to  SPL's  of  165  dB  or  less. 

Earnshaw  (i860)  used  a  piston  as  the  wave  source.   The  piston  then 
represented  a  moving  boundary  condition  and  introduced  into  the  integral 
solution  of  Eq.  (3-l)  a  term  proportional  to  the  piston  displacement. 
Our  solution,  since  it  is  formed  using  a  stationary  boundary,  does  not 
contain  a  source  displacement  term.   In  any  event,  for  propagation 
distances  large  compared  to  the  piston  displacement,  the  piston  displace- 
ment term  contributes  negligibly  to  the  solution. 
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The  source  velocity,  I.e.,  the  velocity  at  x=0,  is  given  by  g(t). 
The  value  of  particle  velocity  at  the  source  for  a  particular  wavelet 
is  then 


u  =  g(t) 


(3-2) 


We  will  solve  Eq.  (3-l)  subject  to  the  condition  that  u  remain 
constant.  Identifying  the  parameter  cp  as  the  time  that  the  wavelet 
of  particle  velocity  u  left  the  origin,  we  may  write 


u  =  g(cp) 


and 


I'  ~  ■  /' 

Jo  yep 


dx*  =   /   (c  +pu)dt" 

'cp 


(3-3) 


(3-M 


Substituting  Eq.  (3-3)  into  Eq.  (3-^)  and  performing  the  integration, 
we  obtain  the  result 


x  =  [cQ  +  pg(cp)](t-cp) 


(3-5) 


Equations  (3-3)  and  (3-5)  a^e  the  Earnshaw  solution.  The  Earnshaw 
solution  is  more  useful  if  Eq.  (3-5)  is  recast  in  an  approximate  form. 
Substituting  the  retarded  time  t'=t-x/c  into  Eq.  (3-5)>  we  obtain 


t'  =  cp  +  — 
c 
o 


1  +  M£l 


-1 


c 


(3-6) 


2,'q,  .  £MM 

c 
o 


3C 


In  Eq.  (3-6)  terms  in  u/c  of  order  greater  than  one  have  been 


ignored.   If  "we  define 


e  =  u  /c     ,  (3-7) 

o'  o  ' 


■where 


u  =  the  peak  value  of  u, 

o  ' 

we  see  that  the  approximation  used  to  obtain  Eq.  (3-6)  is  valid  for 


6  «  1 

For  e-0.1  in  a  sinusoidal  wave,  the  sound  pressure  level  is  about 
175  dB. 

Equation  (3-6)  describes  the  progressive  distortion  of  the 
original  time  base  of  a  signal,  where  the  original  time  base  is 
represented  by  the  parameter  cp  and  the  new  time  base  by  t'.   The 
factor  Pu/c   may  be  thought  of  as  a  time  shift  per  unit  of  propaga- 
tion distance. 

In  summary,  consider  a  source  particle  velocity  relation 
of  the  form  of  Eq.  (3-2).   The  Earnshaw  solution  is 

u  =  g(cp)    ,  (3-3) 

where 

»lt't  £2£M    .  (3.6) 

c 
o 

Equation  (3-3)  may  be  used  to  find  the  value  of  u  at  the  field  point  (x,t) 
provided  Eq.  (3-6)  can  be  solved  for 

cp  =  cp(x,f  )    .  (3-8) 
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The  requirement  to  find  an  explicit  solution  for  cp  in  terms  of  x  and  t' 
is  often  difficult  to  meet  -when  the  Earnshaw  solution  is  used  analytically. 
In  this  dissertation,  we  shall  use  a  computer  algorithm  based  in  part  on 
Eq.  (3-6).  The  numerical  evaluation  of  Eq.  (3-6)  provides  a  calculation 
of  the  distortion  of  the  signal  time  base  and  avoids  the  requirement 
posed  by  Eq.  (3-8). 

b.   Shock  Propagation  Relations 

In  most  instances,  use  of  Eq.  (3-6)  leads  to  multivalued 
waveforms.  The  onset  of  mult i value dness  indicates  the  formation  of 
a  shock.  At  the  point  of  shock  formation  the  Earnshaw  solution  is  no 
longer  valid.  Pressure,  density,  particle  velocity,  enthalpy,*  and 
entropy  on  either  side  of  the  shock  and  the  propagation  velocity  of  the 
shock  are  described  by  use  of  the  Rankine-Hugoniot  shock  relations  Lsee, 
for  example,  Thompson  (1972)J.  Before  stating  these  relations,  -we  will 
establish  the  geometry  and  variables  that  will  be  used  in  the  ensuing 
discussion.  Consider  a  plane  shock  wave  propagating  in  a  perfect  gas. 
The  shock  velocity  is  designated  V;  variables  downstream  (in  front  of) 
the  shock  are  identified  by  the  subscript  1  and  those  upstream  (behind) 
the  shock  by  the  subscript  2.  The  variables  that  we  will  use  are 
summarized  in  Fig.  (3-l)« 


Specific  enthalpy  is  defined  as 


h  =  e  +  Pt 


where  e  is  the  specific  internal  energy,  P  the  absolute  pressure, 
and  t  the  specific  volume  (see,  for  example,  Thompson  (1972)). 
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w2—  V  -  u 

Shock 


Figure  3-1 

A  summary  of  the  variables  used  in  the  Rankine-Hugoniot  relations 

w  is  the  relative  shock  velocity  and  r  is  the  specific  volume 

All  other  variables  are  as  previously  defined 


In  general,  the  wave  flow  fields  on  either  side  of  a  shock  are  varying 
fields.  Unless  stated  otherwise,  in  the  applications  of  the  Rankine- 
Hugoniot  relations  that  we  shall  consider  the  wave  fields  on  either  side 
of  the  shock(s)  will  be  varying.   The  Rankine-Hugoniot  relations  may 
be  stated  as  follows : 


the  continuity  relation 


P2W2  "  P1W1  =  °    ' 


(3-9) 


the  momentum  relation 


the  entropy  relation 


-  P1  +  p2w2     -  plvl     =  0  ,     (3-10) 


1/2        2\ 
the  energy  relation  h2  ~  hl  +  2   \W2  ~W1   /  =  °        »    and  (3-11) 


S£    -    8       >    0 


(3-12) 
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Our  first  task  is  to  calculate  the  magnitude  of  the  entropy 
change  at  the  shock.  The  first  step  in  this  task  consists  of  combining 
Eqs.  (3-9))    (3-10)*  and  (3-ll)  into  a  single  relation  that  contains  only 
thermodynamic  variables.  This  relation,  obtained  after  some  manipulation, 
is 

h2  -  h1  =  t^-p^  + 1  (vV^VV     •  (3"15) 

The  enthalpy  and  specific  volume  upstream  of  the  shock  may  be  described 
in  terms  of  the  downstream  enthalpy  and  specific  volume  by  expanding 
h(P, s)  and  t(P,  s)  in  Taylor  series,  for  example, 

2 


H-\*Ww*$\te^-*-- 


...♦£  (vp! 


,♦4}  ML...       oiu) 


■■■*^p)    (vsi)(p2-V  *  •••      • 

The  series  expressions  may  be  evaluated  explicitly  using  the  perfect 
gas  relations  and  the  identities 


*  -  I  J  (5-15) 

P 


and 


t  -  |  )    .  (3-16) 


3^ 


Substituting  the  evaluated  series  into  Eq.  (3-13)  and  solving  the 
resulting  expression  for  (sp-s  ),  one  obtains 

s2  -  a     =  £*§-  +  o(5^  +  ...    ,  (3-17) 

67 

where 

R  =  the  gas  constant. 
Hence  the  change  in  entropy  across  the  shock  is  of  the  order  of  the 
shock  strength  cubed. 

Until  now,  -we  have  assumed  simple  wave  flow.   Strictly 
speaking,  the  formation  of  a  shock  invalidates  this  assumption.   The 
shock  represents  a  discontinuity  in  the  wave  flow  field.   In  general, 
shocks  are  subsonic  with  respect  to  the  flow  behind  the  shock*  and, 
hence,  a  portion  of  the  wavelets  overtaking  the  shock  is  reflected 
from  it.  The  presence  of  the  reflected  wave  negates  the  simple  wave 
assumption.   However,  if  the  reflected  wave  is  very  weak,  we  may 
ignore  it  and  thereby  preserve  the  simple  wave  assumption.  A  measure  of 
the  strength  of  the  reflected  wave  may  be  gained  by  investigating  the 
magnitude  of  the  discontinuity  in  the  Riemann  invariant  s  at  the  shock. 
Lighthill  (1956)  has  made  such  an  investigation  and  reports  that  the 
change  in  the  quantity  s  is  of  the  order  of  the  change  in  entropy. 
We  see  from  Eq.  (3-17)  that  the  change  in  entropy  is  very  small  for 


That  this  is  true  for  weak  shocks  will  be  apparent  from  the  subsequent 
examples  contained  here.   The  validity  of  the  comment  as  applied  to 
stronger  shocks  is  substantiated,  for  example,  by  Thompson  (1972). 
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weak  shocks  and  can  safely  be  ignored;  similarly,  the  reflected  wave 
can  also  "be  ignored.  Therefore,  in  this  work  simple  wave  flow  will  be 
taken  to  exist  on  either  side  of  the  shock. 

Having  established  the  magnitude  of  the  entropy  change 
associated  with  shock  propagation,  we  desire  an  expression  for  the 
shock  propagation  velocity.   Such  a  relation  may  be  obtained  from  the 
Rankine-Hugoniot  relations.  Equations  (3-9)>  (3-10)>  an(i  (3-H)  may 
be  combined  to  yield  the  relation 

v2  -  pvCug-x^)  -  cx2  =  0    .  (3-18) 

When  Eq.  (3-l8)  is  solved  for  the  shock  speed  V,  we  obtain 


:i  +  Ui  +  f(vV+§!(!^L+... 


— — -  (3-19) 


Substituting  Eq.  (2-3)  evaluated  in  the  downstream  region  into  Eq.  (3-19)> 
we  obtain  the  relation 


T.viW^^(-¥j) 


-l 

♦  ...  •  (3-20) 


Equation  (3-20)  gives  the  propagation  velocity  of  the  shock.  Recall 
that  in  Eq.  (3-6)  terms  of  order  greater  than  e  were  ignored.  When 
we  apply  this  criterion  to  Eq.  (3-20),  we  obtain 


/  shock 


36 


Equation  (3-2l)  is  the  -well-known  -weak-shock  result  that  the  shock 
velocity  is  approximately  the  mean  of  the  propagation  velocities 
immediately  upstream  and  downstream  [see,  for  example,  Thompson  (1972)] . 

Equation  (3-2l)  has  units  of  velocity  and,  in  order  to  be 
compatible  with  Eq.  (3-6),  we  need  to  recast  it  in  units  of  time  shift 
per  unit  distance.  This  recasting  is  easily  accomplished  by  inverting 
Eq.  (3-21),  expanding  the  right  hand  side  in  powers  of  (u  +up)/c  and, 
as  before,  discarding  terms  above  the  first  order.   In  terms  of  the 
retarded  time  t '  we  then  have 


sr)„      =-t4[V^       •         (5-22) 

/shock     2c 
o 


The  quantities  up  and  u..  are  obtained  from  the  Earnshaw  solution 
results  immediately  in  front  of  and  behind  the  shock.   The  values  of 
u„  and  u  used  in  Eq.  (3-22)  provide  the  means  of  connecting  the 
Earnshaw  solution  on  one  side  of  the  shock  with  the  Earnshaw  solution 
on  the  other  side. 

c.   Application  of  Weak-Shock  Theory 

In  Subsections  (a)  and  (b)  we  have  derived  the  basic  mathematical 
tools  used  in  weak-shock  theory.  These  tools  are  the  Earnshaw  solution 
and  the  shock  arrival  equation.  We  must  also  provide  a  source  relation 
of  the  form  of  Eq.  (3-2).  Equation  (3-17)  also  is  used  in  weak-shock 
calculations,  but  in  an  indirect  way.   Equation  (3-17)  shows  that  there 
is  dissipation  at  the  shock.   The  loss  mechanism  is  built  into  Eq.  (3-2l), 
however,  and  we  shall  never  have  to  use  Eq.  (3-17)  explicitly. 
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Rudnick  (1952)  made  explicit  use  of  Eq.  (j3-17)  and  for  the  case  of 
sawtooth  -wave  propagation  obtained  the  -weak-shock  theory  result. 

We  shall  apply  weak-shock  theory  to  two  problems .  In  the 
first  problem,  only  the  Earnshaw  solution  portion  of  the  theory  is 
used.  In  the  second  problem  the  full  weak-shock  theory  is  demonstrated, 
l)   The  First  Example  Problem* 

Consider  a  sound  wave  whose  shape  at  its  source  is  as  shown 
in  Fig.  (3-2). 


x  =  O,  t'=  t  s  <P 


Figure  3-2 
The  source  waveform  for  the  first  example- 


The  source  waveform  is   described  analytically  by  the  relations 


u  = 


u 

-2.  t 

T  t 

o 


=  u 


0  ^  t  ^  T 


T  S  t 
o 


(3-23) 


-if 


Rayleigh  (1910)  considered  an  example  similar  to  our  first  example 
problem.   Rayleigh' s  example  was  treated  as  an  initial  value  problem. 
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Equation  (3-23)  is  the  source  relation  [see  Eq.  (3-2)]  for  this  example. 
Substituting  Eq.  (3-23)  evaluated  at  t=cp  into  Eq.  (3-6),  "we  obtain 


=  t* 


1- 


Pu  x 
o 


1-1 


T  c 
o  o 


for 


and 


0  ^  cp  ^  T     , 
o 


pu  X 


cp  =  tf  + 


(3-24) 


(3-25) 


for 


cp  ^  T 


Equations  (3-24)  and  (3-25)  are  the  explicit  solution  of  cp  in  terms  of 
x  and  t'  indicated  earlier  in  general  terms  by  use  of  Eq.  (3-8) •   The 
solution  for  u  as  a  function  of  cp  -where  cp=cp(x,  t1)  is  then 


for 


u(x,f)  =  ^f 

o 


1  - 


Pu  X 
Q 


1-1 


2 
T  c 
o  o 


u  ex 

0  ^  t*  ^  T  -  -^5-   , 


(3-26) 


and 


u(x,t')  =  u 


(3-27) 


for 


t'  S  T  - 


u  px 
o 


°     c2 

o 


Use  of  Eqs.  (3-26)  and  (3-27)  describes  the  deformation  of  the  wave 
shown  in  Fig.  (3-2)  for  some  propagation  distance  x.  Note  that  at 
x=0  Eqs.  (3-26)  and  (3-27)  reduce  to  Eq.  (3-23). 

The  distortion  evaluated  by  use  of  Eqs.  (3-26)  and  (3-27) 
is  shown  in  Fig.  (3-3)- 


Uo 

X  :0 

A 

T0                                  t 

(a) 


0  <  x  <x 


(b) 


(c) 
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X    =   X 


t' 


Figure  3-3 

The  progressive  distortion  of  the  Fig.  3-2  waveform 

(a)    x  r  0,      'b)    0<  x<  7,      (c)    x  =  x" 


Figure  (3~3a)  repeats  the  source  waveform  shown  in  Fig.  (3-2).   In 
Fig.  (3-3"b)  the  wave  has  been  distorted.   Note  from  this  example  that 
straight  line  sections  remain  straight  as  they  are  distorted;  this  is  a 
general  result.   Figure  (3-3c)  shows  the  waveform  at  the  shock  formation 
distance  x.   The  quantity  x  may  be  calculated  using  Eq.  (3-26)  by  noting 
that  u  becomes  discontinuous  at  the  point 


x 


c  T 
o  o 


Pu 


(3-28) 


Equation  (3-28)  is  also  a  general  result  applicable  to  straight  line 
sections  not  immediately  preceded  or  followed  by  a  shock. 


If  we  continue  to  apply  the  Earnshaw  solution  to  the  above 
problem  for  x  >  x  we  would  obtain  multiple  valued  solutions.   Hence, 
further  analysis  requires  use  of  the  weak-shock  relations.  In  the  next 
example  problem  the  use  of  the  weak-shock  relations  will  be  demonstrated 
for  a  case  containing  two  shocks . 

2)   The  Second  Example  Problem 

In  Fig.  (3-4a)  we  show  the  source  waveform  for  the  second 
example  problem.   Note  that  the  slopes  of  the  segments  AB  and  CD  are 
equal.   In  accordance  with  results  obtained  by  use  of  Eq.  (3-28)  both 
segments  will  form  shocks  at  the  same  value  of  x.   Having  equal  slopes 
for  the  segments  AB  and  CD  is  merely  a  convenience  to  facilitate  the 
subsequent  calculations.   Had  the  initial  slopes  of  the  segments  AB  and 
CD  been  different,  the  calculations  would  have  proceeded  in  a  similar, 
albeit  slightly  more  complex,  fashion. 


uD(o) 


UD(X) 


un(x) 


t;(*)  t^(x)  t'pix)        t'E(x) 


t' 


Figure  3-4 

The  second  example  problem  waveforms 

(a)    x  =  0,      (b)    0<x  <x,      (c)    x  =  x 


in 


In  Fig.  (3-^ti)  the  wave  has  been  distorted  hut  shocks  have  not  yet 
formed.  Figure  (3-1+c)  shows  the  -waveform  at  the  point  of  shock  formation, 
Henceforth,  -we  shall  refer  to  the  segments  AB  and  CD  as  the  first  and 
second  shocks,  respectively.  The  notation  that  we  shall  use  is  as 
follows : 

a)  The  capital  letters  A,  B,  C,  D,  and  E  indicate  points  of  slope 
change  of  the  waveform  for  some  propagation  distance  x.  The 
point  F  is  the  point  of  intersection  of  the  segment  BC  and 
the  t '  axis . 

b)  The  coordinates  of  a  particular  point,  such  as  point  A,  for 
some  value  x  are  (u.(x),t'(x)). 

c)  Values  of  u  along  the  segments  BC  and  DE  for  some  value  x 
are  indicated  by  the  notation  il  (x,tT)  and  u  (x,f). 
The  distortion  of  the  waveform  shown  in  Fig.  (3-ka.)   to  the 

point  of  shock  formation  proceeds  exactly  as  indicated  for  the  first 
example  problem.   Let  us  assume  that  the  wave  has  propagated  to  the 
point  x  as  shown  in  Fig.  (3-3c)  and  that  we  wish  to  compute  the  subse- 
quent waveform  distortion.  For  this  calculation  we  shall  take  the  space 
origin  to  be  the  point  at  which  the  shocks  form.  The  time  of  arrival  of 
a  particular  wavelet  at  the  origin  is  indicated  by  the  parameter  cp.  At 
the  origin  cp,  t,  and  t'  are  equal.   The  source  relations  for  this  problem 
are 


uBC(o,t')=  M^Cyfo)  -  t'J  (3-29) 


for 


tA'(o)  g  t1  g  tc'(o) 


kz 


and 


for 


where 


UpgCo^')  =  %EttE'(o)  -  t»]  (3-30) 


tc'(o)  §  f  =s  tE'(o) 


^BC  =  V(0)-tA'(0)     >  (5"51) 


and 


^E  =  v(o)  -  t  '(o)    '  (5-32) 


E  v  '  C 

Substituting  in  turn  Eqs .  (3-29)  and  (3-30)  evaluated  at  t=cp  into 
Eq.  (3-6),  we  obtain  expressions  for  cp  as  a  function  of  (x, t').   The 
solutions  for  u  as  a  function  of  cp(x,  t')  are  then 


M^Ctp'Co)  -  f  ] 


V^'>  =  ^  %Vx  C-5?) 


1  + 


2 
c 
o 


for 


tA'(x)  §  t»  ^  tc'(x) 


^ 


and 


VY(0)it,] 

)E3 


V'^''-       PMdfx  <^> 


i  +  — r 

c 
o 

for 

tc'(x)  *  f  £  tE'(x) 

At  this  point  we  do  not  know  the  value  of  t  '(x)  and  t  '(x)  because 
these  quantities  must  he  calculated  using  Eq.  (3-22).  The  quantity 
t  '(x)  remains  a  constant  equal  to  t  ' (o)  because  the  particle  velocity 

ill  ill 

at  point  E   is   zero. 

We  are  now  ready  to  compute  the  motion  of  the  shock  fronts. 

This  computation  is  performed  using  Eq.  (3-22).   In  this  problem  we 

have  identified  the  times  of  occurrence  of  the  shocks  as  tA'(x)  and  t  '(x) 

Hence,  in  using  Eq.  (3-22),  we  adopt  the  notation 


dx      dx 

shock  1 


(3-35) 


and 


(3-36) 


dtc'(x)   dV 

^     "  ^  'shock  2 

where  shock  1, 2  denotes  the  first  or  second  shock,  respectively. 

Evaluating  Eq.  (3-22)  at  the  first  shock,  we  obtain 


dt  '(x) 

2c 
O 


kk 


At  the  first  shock  the  velocity  downstream  of  the  shock  is  zero;  hence, 
the  quantity  u  does  not  appear  in  Eq.  (3-37) •  Substituting  Eq.  (3-33) 
evaluated  at  t  '(x)  into  Eq.  (3-37)>  ™e  obtain 


dtA'(x) 
dx 


[t  '(0)  -  V(x)] 


2c 


"^C 


A 


1  + 


PMBC 


(3-38) 


The  solution  of  Eq.  (3-38)  subject  to  the  boundary  condition 
t  '(x)=t  '(0)  at  x=0  is 


tA'(x)   =   tA'(0) 


1  + 


^ 


x 


11/2 


+  y(o)(i 


i  + 


PMBC 


11/2' 


o 


(3-39) 


The  motion  of  the  second  shock  is  somewhat  more  complicated 
because  the  velocity  downstream  of  this  shock  is  not  zero.   Substituting 
Eqs.  (3-33)  and  (3-3^)  evaluated  at  t  '(x)  into  Eq.  (3-22),  we  obtain 


dx 


2c 


W(0)  -  V<x»    •WV*05  -  Vix)) 


i  + 


PMK 


1      + 


^MDEX 


(3-^0) 


^ 


The  solution  of  Eq.  (3-^°)  subject  to  the  boundary  condition  t  '(x)=t  ' (0) 

L>        C 

at  x=0  is 


tc'(x)  =  tc'(0) 


1  +  "^  "bc* 

C 
O 


1/2 


1  + 


PMDI 


,x 


.11/2 


tF'(° 


1  + 


PMDI 


/^ 


?75 


o 


1/2 


+    1  + 


PMDI 


,X 


jm 


(3-^1) 


tE'(0)pMDEX 


1  + 


PMBC 


W 


1/2 


^1 


1  + 


PMBCX 


+  1  + 


PMDI 


Equations  (3-39)  and  (3-^1)  give  the  retarded  times  of  occurrence  of 
the  two  shocks  for  any  propagation  distance  x  up  to  the  point  of  shock 
merger. 

Having  solved  for  tfl'(x)  and  t  '  (x)  "we  may  now  calculate  the 
values  of  particle  velocity  along  the  segments  BC  and  DE.   Two  values 
of  particular  interest  are 


ub(x)  =  v(x,V(x))     ' 


(3-te) 


and 


VX)  =  UDE(X^C,(X)) 


(3-^3) 


h6 


The  expression  for  U-n(x)  turns  out  to  be  particularly  simple,  i.e. 


Vx) 


V°) 


1  + 


PMBCX 


W2 


2 
c 
o 


(3-W) 


The  result  for  U-pXx)  is  similar  but  more  complicated  and  for  the  sake  of 

brevity  -we  shall  not  present  it  here.  The  amplitude  decay  function 

/  .      2\1/2 

ll+piVL  x/c      I         in  Eq.    (3-h-k)   is  of  the  same  form  as  that  for  N  waves 

[see,  for  example,  Blackstock  (1972)],  a  not  surprising  result,  since 

the  triangular  section  ABCF  of  the  original  -waveform  (Fig.  3-^c)  is 

really  half  of  an  N  -wave. 

Let  us  now  consider  the  merger  of  the  two  shock  waves.   In 

this  problem  it  turns  out  that  for  some  value  of  x  we  obtain  the  result 

tA'(x)  =  tc'(x)    .  (5-^5) 

The  value  of  x  satisfying  Eq.  (3-^-5)  can  be  found  analytically  or  in  a 
specific  case  by  numerical  evaluation.   Physically,  the  equality  of  the 
shock  arrival  times  indicates  that  the  two  shocks  have  merged.   In  this 
example,  the  merger  results  when  the  second  shock  overtakes  the  first. 
Following  the  shock  merger,  the  resultant  waveform  appears  as  in 
Fig.  (3-5d.).   The  equations  describing  the  retarded  time  of  arrival  of 
the  merged  shock  and  the  value  of  particle  velocity  behind  the  merged 
shock  are  similar  to  Eqs.  (3-39)  and-  (3-^-*0>  respectively. 

A  computer  program  was  written  to  evaluate  the  results  for 
this  problem.   The  waveform  depicted  in  Fig.  (3_1+c)  was  used  as  the 
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source  -wave.   The  computations  were  done  in  two  stages.  First  Eqs .  (3-39)> 
(3-kl)t    (~5-Kh),  a.nd.   the  relation  for  u^(x)  not  given  here  were  used  up 
to  the  point  of  shock  merger.  Following  shock  merger,  equations  similar 
to  Eqs.  (3-39)  an(i  (3-^*0  "were  used.   The  results  appear  in  Fig.  (31-!?)* 
The  distances  indicated  on  each  drawing  are  distances  in  feet  with 
respect  to  the  original  shock  formation  point  for  this  problem.  The 
values  of  the  remaining  parameters  at  the  original  shock  formation  point 
are  t  '=1.0  msec,  t  '=1.25  msec,  t  '=2.0  msec,  u=5  ft/sec,  and 
u=15  ft/sec. 
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Figure  3-5 
The  second  example  problem  (cont'd) 
The  source  waveform  in  the  shifted  coordinate  system. 
The  source  waveform  after  propagating  25  ft. 
Just  prior  to  shock  merger. 
The  simplified  waveform  after  shock  merger. 
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The  two  example  problems  illustrate  several  -well-known  effects 
that  take  place  in  finite -amplitude  waves.  These  effects  are  as  follows: 

a)  Nonlinear  effects  cause  waves  to  distort,  i.e.,  positive 
sloping  sections  form  shocks,  and  the  slopes  of  negatively 
sloped  sections  gradually  decrease. 

b)  Shocks  exhibit  relative  velocity  with  respect  to  one  another. 
Hence  shocks  may  merge  as  the  wave  propagates . 

c)  The  wave  is  attenuated  even  though  attenuation  is  not 
explicitly  included  in  Eq.  (3-2l).   The  reason  for  this,  as 
discussed  earlier,  is  that  attenuation  is  "built  into" 

Eq.  (3-2l).   Physically,  the  attenuation  is  a  consequence  of 
the  irreversible  nature  of  shock  propagation.   This  point  was 
noted  in  the  discussion  following  the  derivation  of  Eq.  (3-1?) • 

d)  The  shape  of  the  wave  is  usually  significantly  altered  as  it 
propagates.   Hence,  we  expect  the  frequency  spectrum  to  be 

a  function  of  propagation  distance. 

e)  Straight  line  segments  distort  in  such  a  fashion  that  they 
remain  straight.  The  slope  M  is  related  to  the  initial  slope 
M  by  an  equation  of  the  general  form 

M 

1 °_ 

1     2 
c 
o 

Positively  sloping  sections  form  shocks.   The  shock  formation 

distance  for  sections  not  immediately  preceded  or  followed  by 

a  shock  in  general  is 

2 
x  =  —  (5-^7) 
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3.   A  General  Technique  for  Modeling  Arbitrary  Waveforms  "by  Straight 
Line  Sections 

The  method  used  in  the  preceding  example  may  be  generalized  for 

use  with  any  sound  wave  modeled  by  straight  line  sections.   In  this 

section,  the  general  solution  for  a  plane  shock  wave  preceded  and  followed 

by  a  straight  line  velocity  field  of  arbitrary  slope  is  calculated.  The 

waveform  at  the  space  origin  is  shown  in  Fig.  (3-6). 
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Figure  3-6 
A  plane  shock  wave  propagating  in  a  disturbed  medium 


The  notation  that  we  use  in  this  section  follows  the  same 
rules  as  the  notation  used  earlier  in  the  second  example  problem  in 
Section  2  of  this  chapter.   Expressions  for  the  retarded  time  of  shock 


arrival  t  ' (x)  and  for  the  particle  velocities  in  front  of  and  behind 
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the   shock  are  obtained  exactly  as  they  -were  obtained  in  Section  2. 
These  expressions   are 


tB'(x)   =   tB'(0)yi   -  A^x/l^    A 


CD 


ACD   "  V 


/l   -  A^l   -  ACDx   -1  ♦  A^x] 
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W'0'  -p°o   V°'r 


ACD   "  AAB 


x/l   -  k„xJl   -  A„x  -1  +  A„  xl 

v  AB    v  CD  CD   J 


V*>  = 
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and 


2.-1 


uc(x)   = 
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1   "  ACDX 
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■where  A       and  A-^  (closely  related  to  M.      and  WL   )   are  defined  by 
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/  vo)  -  vo) 

>  '  \tA'(o)  -  tB'(o) 


(3-51) 


and 


p     /  uc(0)  -  V0) ' 

A0D  "  ~  ■    (^(0)    -  tD'(0) 

O  X 


(3-52) 
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Equations  {j-kQ)   through  (3-52)  are  applicable  until  the  propagating 

shock  coincides  with  one  of  the  end  points  of  the  -wave  section  shown 

in  Fig.  (3-6).  Additionally,  these  end  points  could  be  overtaken  by 

some  other  shock  outside  the  region  defined  by  Fig.  (3-6).  Points 

labeled  A  and  D  in  Fig.  (3-6)  "will  themselves  be  shifted  as  described 

by  Eq.  (3-6).  The  earliest  time  of  coincidence  of  t  '(x)  "with  either 

B 

t.'(x)  or  t  '  (x)  can  be  calculated.  Coincidence  -with  other  shocks 

propagating  from  outside  the  region  shown  in  Fig.  (3-6)  can  also  be 

computed.  Once  the  earliest  coincidence  distance  is  reached  a  new  set 

of  conditions  could  be  established  and  the  solution  continued.  While 

the  method  appears  to  be  extremely  unwieldy,  it  probably  could  be  made 

to  -work,  -within  the  limitation  that  straight  line  modeling  be  used. 

Such  an  attempt,  in  a  very  simple  case,  has  been  reported  by  Hawkings  (1970) 

and  was  discussed  in  Chapter  2.  Hawkings'  computer  method  produces 

results  that  exhibit  shock  propagation  characteristics  similar  to  those 

shown  in  Fig.  (3-5). 

The  complexity  of  solutions  of  the  type  typified  by  Eqs.  (3-h8) 
to  (3-50)  and  the  great  difficulty  in  adapting  them  to  individual  cases 
convinced  us  that  a  numerical  approach  to  the  problem  would  be  easier. 
This  approach  is  described  in  Chapter  k. 

In  Chapter  2  we  mentioned  that  use  of  weak-shock  theory 
produced  results  equivalent  to  some  of  the  results  of  Burgers  (1972). 
We  will  now  demonstrate  this  fact.   In  the  ensuing  discussion  let  us 
assume  that  the  propagation  distance  x  is  very  large.   Noise  waveforms 
at  very  large  values  of  x  will  have  deformed  into  a  series  of  shock 
fronts  separated  by  negatively  sloped  (c\i/oV<D)  straight  line  sections. 
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The  noise  field  -will  then  appear  as  a  train  of  randomly  spaced  shock 
waves.   The  velocity  at  any  point  on  one  of  the  negatively  sloped 
sections  will  be  described  by  an  equation  of  the  form  of  Eq.  (3-I49). 
That  is,  the  velocity  on  a  negatively  sloped  section  is  of  the  form 


u(x,t»)  =  2 .       (3.53) 

1  "  AABX 


where  the  constants  in  Eq.  (3-53)  are  as  defined  for  Eq.  (3-51).  The 
slope  of  Eq.  (3-53)  for  constant  x  is 


For  very  intense  waves  at  long  distances  from  the  source  A  x|»l,  and 
we  obtain 

|r— 7^    •  (3-55) 

/  X 


This  result  which  is  equivalent  to  results  reported  by  Burgers  (1972) 
shows  that  slopes  of  the  continuous  portions  of  a  finite -amplitude  wave 
under  the  aforementioned  conditions  are  all  equal  and  are  proportional 
to  x   .  Moreover,  the  slopes  are  independent  of  the  original  slopes. 
The  asymptotic  form  of  a  noise  wave  is  then  a  series  of  randomly  spaced 
shocks  separated  by  straight  line  sections  of  equal  slope.   Since  the 
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intershock  sections  are  of  equal  slope,  the  shock  velocities,  though 
individually  varying,  remain  constant  until  merger  with  another  shock. 

k.        Impedance  Relations 

Both  Eqs.  (3-6)  and  (3-2l)  are  in  terms  of  particle  velocity. 
However,  we  actually  measure  pressure.   In  the  continuous  portion  of 
the  wave,  pressure  can  be  easily  related  to  particle  velocity.  Making 
use  of  Eq.  (2-2)  and  the  isentropic  perfect  gas  relations,  we  obtain 
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27 

r 


(3-56) 


where  the  subscripts  1  and  2  are  used  to  indicate  the  beginning  and 
ending  values  of  the  variables  in  an  isentropic  process.   Substituting 
Eq.  (2-3)  into  Eq.  (3-56),  we  obtain  the  result  for  outgoing  waves, 
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Expanding  Eq.  (3-57)  in  powers  of  (up-u  )/c  yields  the  result 
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2    1 


P1c1(u2-u1) 
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.  (3-58) 


A  result  similar  to  Eq.  (3-58)  may  be  derived  for  the  pressure 
change  across  a  shock.   Combining  Eqs.  (3-9),  (3-10),  and  the 
relative  shock  velocity  definition  (see  Fig.  (3-l))>  we  obtain 


U2'U1  =  [(P2'P1)(T2"T1) 


1/2 


(3-59) 


The  specific  volume  t(s,P)  is  expanded  in  a  Taylor  series  and  this 
series  is  substituted  into  Eq.  (3-59).   The  binomial  theorem  is  then 
used  to  expand  the  result,  yielding 
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By  use  of  the  reversion  relations  [see,  for  example, Chemical  Rubber  Company 
Handbook (196^)]  Eq.  (3-6O)  may  be  written 


P2   '  Pl  =  P1C1(VV 


p  <vV    p2  (u2-ui)2 


2 


'1 


■     (3-61) 


Hence,  by  comparing  Eqs.  (3-58)  and  (3-6l)  we  may  make  the  following 
statement  concerning  a  simple  wave  transition  and  a  shock  wave  transition, 
The  series  solutions  relating  the  change  in  pressure  to  the  change  in 
particle  velocity  in  either  transition  agree  in  form  for  terms  through 
the  second  power  in  the  change  in  particle  velocity.  Moreover,  if  the 
magnitudes  of  the  transitions  in  the  two  cases  are  equal,  the  series 
are  equal  through  the  second  term.   The  error  involved  in  discarding 
terms  above  the  first  order  in  the  change  in  particle  velocity  may  be 
easily  estimated  by  evaluating  the  magnitude  of  the  discarded  terms. 
In  general,  for  waves  that  may  be  analyzed  by  weak-shock  theory,  this 
error  is  small  and  noncumulative.   Hence,  for  the  remainder  of  this 
dissertation,  we  shall  use  the  linear  impedance  relation 


P  _  p  =  p  c  u    , 
o    00 


(3-62) 
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where  we  have  replaced  the  2  variables  by  a  general  unsubscripted 

variable  and  the  1  variables  by  their  undisturbed  values  P  ,  p  ,  c  and 

o   o   o 

u  =0. 
o 

5.   Tube  Wall  Effects 

In  Chapter  2  we  mentioned  that  McKittrick  et  al.  (1967)  and 
Pernet  and  Payne  (1969)  had  observed  asymmetric  waves  propagating  in 
plane  wave  tubes.   The  source  waveforms  for  both  experiments  were  sinu- 
soids. Using  weak-shock  theory,  we  would  predict  sawtooth  waves. 

Both  McKittrick  et  al.  and  Pernet  et  al.  concluded  that  the 
asymmetry  was  due  to  dispersion  caused  by  the  tube  wall  boundary  layer. 
A  typical  asymmetric  waveform  reproduced  from  the  paper  by  McKittrick 
et  al.  is  shown  in  Fig.  (3-7)- 


Figure  3-7 
A  typical  shock  wave  in  a  plane  wave  tube 


Coppens  (1971)  predicted  that  asymmetry -will  also  be  observed  before 
shock  formation  and  based  on  theoretical  considerations  concluded 
that  dispersion  is  the  cause. 

Since  our  experimental  conditions  are  very  similar  to  those 
under  which  both  McKittrick  et  al.  and  Pernet  et  al.  worked,  it  is 
important  to  include  both  tube  wall  dispersion  and  attenuation  in  our 
analysis.   In  a  wide  tube*  the  approximate  expressions  for  the  phase 
velocity  c  and  attenuation  coefficient  a   obtained  from  the  Kirchhoff 
theory  are 


CD 


1  + 


7-1 
v/Pr" 


% 


(3-63) 


and 
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(3-64) 


A  "wide  tube"  [the  terminology  is  due  to  Weston  ( 1953)1  is  one  whose 

boundary  layer  thickness  A  is  small  compared  to  the  tube  radius.  We 

shall  use  Weston's  definition  of  boundary  layer  thickness,  the  distance 

from  the  tube  wall  to  the  point  where  the  principal  maximum  in  particle 

l/2 
velocity  occurs,  i.e.,  A  =  2.2838(2v/oo)  '  ,  where  00  is  the  angular 

frequency.  A  further  property  of  a  wide  tube  is  that  it  must  not  be 

so  wide  that  losses  in  the  fluid  mainstream  are  comparable  to  the  wall 

losses.   This  latter  requirement  may  be  stated  quantitatively  by 

comparing  the  expressions  for  the  boundary  layer  attenuation  coefficient 

and  for  the  unbounded  medium  attenuation  coefficient ;  one  obtains 

kR  «  T)/kA,  where  r\   is  a  constant  determined  by  the  environmental 

conditions.  At  72°F  T)  equals  approximately  J>.6.      Hence,  a  wide  tube 

2 
is  one  that  meets  the  combined  requirements  A  «  R  «  rj/k  A. 
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Inspection  of  Eq.  (3-63)  shows  that  the  speed  of  small-signal  -waves 
in  a  tube  is  less  than  c  ;  the  speed  approaches  c  only  at  very  high 
frequencies.  Hence,  any  small-signal  nonsinusoidal  -waveform  will  he 
distorted  as  it  propagates  in  a  pipe  because  the  higher  frequency 
components  travel  faster  than  the  lower  frequency  components.  The 
differences  in  phase  velocity  are  small  but  their  effect  is  cumulative 
with  distance.  Use  of  Eq.  (3-63)  allows  us  to  calculate  the  cumulative 
distortion  due  to  dispersion  for  small-signal  processes.  In  a  nonlinear 
problem  the  distortion  due  to  dispersion  is  less  readily  calculated 
because  growth  of  the  various  harmonics  and  their  intermodulation 
products  is  a  continuous  process. 

Equations  (3-63)  an<^  (3-6*0  may  be  used  to  calculate  approximate 
values  of  c  and  OC   for  noncircular  ducts  if  we  replace  R  by  D/2  where  D 
is  the  hydraulic  diameter  [see,  for  example,  Thompson  (1972)].   Therefore, 
since  Eqs .  (3-63)  and  (3-6*4-)  are  the  only  portions  of  modified  weak-shock 
theory  that  depend  upon  the  pipe  shape,  the  use  of  the  hydraulic  diameter 
concept  allows  us  to  also  apply  the  theory  to  propagation  problems  in 
noncircular  ducts. 

6.   Modified  Weak-Shock  Theory 

In  this  section  we  shall  propose  a  new  theory  in  which  nonlinear 
effects  embodied  in  weak-shock  theory  and  tube  wall  effects  are  combined. 
The  new  theory  will  be  called  modified  weak-shock  theory. 
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For  the  convenience  of  the  reader,  -we  will  restate  the  Earnshaw 
solution  and  the  shock  arrival  equation.  The  Earnshaw  solution  is 

u  =  g(cp)    ,  (3-3) 

■where 

9  i  f  +  6£i|^    .  (5-6) 

c 
o 

The  shock  arrival  equation  is 

'  o 

Given  a  source  -waveform,  use  of  the  Earnshaw  solution  and  of  Eq.  (3-22) 
allows  us  to  calculate  the  progressive  distortion  of  the  continuous 
portion  of  the  waveform  and  the  movement  of  the  shock  front,  respectively. 
As  we  have  indicated  previously,  our  solution  technique  involves  a 
computer  algorithm  "based  upon  Eqs.  (3-6)  and  (3-22).  The  method  whereby 
we  account  for  tube  wall  attenuation  and  dispersion  is  very  straight- 
forward. We  theorize  that  the  small-signal  tube-wall  attenuation  and 
dispersion  expressions  are  accurate  for  short  distances.  Accordingly, 
we  expect  each  frequency  present  to  be  attenuated  and  dispersed  as 
described  by  use  of  Eqs.  (3-63)  and  (3-64).   The  attenuation  and 
dispersion  relations  may  be  combined  in  the  more  compact  form  given  by 

-a  (l+j)Ax  (3-65) 

C.  '  =  C.  e 
1     1 
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where 

C .   =  the  complex  amplitude  of  the  ith  frequency  component  of 
1 

the  wave, 

C.'   =  the  complex  corrected  amplitude, 

a         =  the  attenuation  coefficient  evaluated  at  the  angular 
i 

frequency  ox  ,  and 
Ax  -   a  small  distance  increment. 
In  the  use  of  modified  weak-shock  theory,  the  transformations  between 
the  frequency  and  time  domain  are  performed  "by  use  of  an  FFT  (fast 
Fourier  transform)  computer  routine.   In  this  routine  the  various 
components  C.  are  assumed  to  be  harmonically  related.   Therefore,  a.  is 
correctly  replaced  by  vn  Q!  ,  where  OC     is  the  coefficient  for  the  funda- 
mental frequency.   Our  computer  algorithm,  which  is  described  in  detail 
in  Chapter  h,    directly  incorporates  Eqs.  (3-6),  (3-22),  and  (3-65) 
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FIGURE  3-8 
A  SCHEMATIC  DIAGRAM  OF  MODIFIED  WEAK -SHOCK  THEORY 
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cast  in  appropriate  forms.   The  method  is  summarized  in  Fig.  (3-8). 

Implicit  in  the  model  that  we  have  called  modified  weak-shock 
theory  are  the  following  assumptions: 

a)  Only  plane  waves  are  present. 

b)  The  wave  amplitude  is  appropriate,  i.e.,  not  too  large,  for 
the  use  of  modified  weak-shock  theory. 

c)  The  pipe  is  a  wide  pipe,  i.e., 

A  «  R  «  -4L- 
k  A 

d)  The  pipe  is  circular,  infinitely  long  and  rigid. 

e)  The  enclosed  fluid  medium  is  a  homogeneous  perfect  gas. 


CHAPTER  h         THE  WAVE  PROPAGATION  ALGORITHM 

1.   Introduction 

This  chapter  contains  the  detailed  computer  techniques  and 
procedures  used  to  calculate  the  distortion  of  a  finite -amplitude  plane 
wave  of  arbitrary  form  propagating  in  a  pipe.   The  reader  is  advised 
that  he  need  not  read  Chapter  h   unless  he  is  interested  in  the  actual 
computational  details.  For  those  not  so  interested,  "we  suggest  examina- 
tion of  Figs,  (h-^-) ,  Oj-5) }   an(i  (^-6)  and  reading  of  the  paragraph  just 
prior  to  the  figures .  This  material  gives  a  general  idea  of  the  distor- 
tion of  a  high  intensity  noise  -waveform  and  prepares  the  reader  for  the 
results  presented  in  Chapter  6. 

In  the  opening  section  of  the  chapter,  the  -weak-shock  theory 
described  in  Chapter  3  is  restated  and  the  pertinent  equations  are  cast 
into  finite  difference  form;  these  difference  equations  are  solved, 
thereby  opening  the  -way  to  the  formulation  of  a  computer  routine  that 
implements  weak-shock  theory. 

The  method  used  to  incorporate  attenuation  and  dispersion 
into  the  basic  weak-shock  algorithm  is  then  considered,  and  the  final 
form  of  the  modified  weak-shock  theory  is  presented.  The  next  portion 
of  the  chapter  is  devoted  to  outlining  the  solutions  of  some  minor 
computational  problems  and  indicating  how  we  selected  values  for 
several  important  program  parameters . 

The  chapter  is  concluded  with  a  discussion  of  the  programs 

that  were  used  to  calculate  power  spectral  densities. 
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2.   Discretization  of  the  Weak-Shock  Theory  Equations 

The  basic  equations  of  -weak-shock  theory  -were  derived  in 
Chapter  3«  For  the  convenience  of  the  reader  we  shall  repeat  the 
equations  here.   Consider  a  plane  wave  whose  time  waveform  at  the 
source  is  given  by 

u  =  g(t)   .  (3-2) 

The  particle  velocity  away  from  the  source  is  then  described  by 

u  =  g(cp)    ,  (3-3) 

where 

cp  =  t'+^|i2    .  (3-6) 

c 
o 

Equations  (3-3)  and  (3-6)  constitute  the  Earnshaw  solution.   The 
Earnshaw  solution  is  used  in  weak-shock  theory  to  describe  the  distortion 
of  the  continuous  portions  of  the  propagating  sound  wave .   Once  shocks 
have  formed,  the  retarded  time  of  occurrence  of  a  particular  shock  is 
described  by  the  differential  equation 

^T1)  =   "  -2-2   k(x,t'(x))  +  u2(x,f  (x))]     .   (3-22) 

/  shock     2c 
/  o 

The  values  of  u  and  u  used  in  Eq.  (3-22)  provide  the  means  of  connecting 
the  Earnshaw  solution  on  one  side  of  the  shock  with  the  Earnshaw 
solution  on  the  other  side. 
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It  is  convenient,  at  this  point,  to  extend  our  interpretation 

of  Eq.  (3-6)  so  that  use  of  Eq.  (3-6)  describes  the  time  shift  between 

t'  and  any  previous  retarded  time  rather  than  the  time  shift  between  t1 

and  the  time  of  generation  of  the  -wavelet.   The  extension  is  implemented 

by  replacing  cp  by  t '    ,  t '  by  t '    ,  and  x  by  Ax  -where  Ax  represents 

the  incremental  propagation  distance  over  which  t'  . .  is  distorted  into 

old 

t'    .   Once  -we  have  made  these  substitutions,  we  obtain 

new  ' 

Pu(t*    )Ax 

V  =  f.  . ^ .  (4-1) 

new     old         2 

c 
o 

Equations  (4-l)  and  (3-22)  constitute  a  system  that  we  shall 
use  to  describe  the  distortion  of  a  source  -waveform  as  it  propagates. 
Our  goal  is  to  determine  the  output  of  this  system  for  an  arbitrary 
input  waveform.   For  simple  cases,  this  determination  can  be  made  by 
analytic  methods,  but,  as  Sections  2  and  3  of  Chapter  3  suggest,  the 
analytic  methods  become  increasingly  cumbersome  for  complicated  waveforms 
Hence,  we  desire  a  computer  method  that  describes  the  continuous  system 
represented  by  Eqs .  (3-22)  and  (4-1).   In  practice  a  continuous  system 
can  be  approximated  by  replacing  the  continuous  equations  by  difference 
equations.   This  is  required  if  -we  wish  to  use  a  digital  computer  to 
calculate  the  system  response  to  various  inputs. 

In  the  system  under  consideration,  the  governing  relations  are 
very  simple  and  the  difference  equations  can  be  written  by  inspection. 
For  higher-order  systems,  the  transition  to  a  difference  format  is  not 
so  simple  and  a  formal  mathematical  technique  is  required  to  ensure  the 
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generation  of  the  proper  equations .   The  state  variable  techniques  used 
by  Ogata  (1967)  constitute  such  a  method.   Proceeding  either  by  the 
state  variables  method  or  simply  by  inspection,  -we  obtain  the  difference 
forms  of  Eqs.  (4-1)  and  (3-22).   These  forms  are,  respectively, 


t'[(k+l)h]  =  t*[kh]  -  pc  "2u[t'(kh)]h  (4-2) 


and 


phc  "' 

t  '[(k+l)h]  =  t  *  [kh] £— 

s  s         d 

where 


u2 


[kh,t  '(kh)]  +  u.[kh,t  '  (kh)]|  ,( 4-3) 


h  =  an  incremental  distance, 

t'[(k+l)h]  =  the  value  of  t'  after  k+1  incremental  steps  of  size  h, 

t  '  =  the  value  of  t'  associated  -with  a  particular  shock, 
s 

t  '[(k+l)h]  =  the  value  of  t  '  after  k+1  steps  of  size  h,  and 
s  s 

u   [(kh,t  '(kh)]  =  the  value  of  particle  velocity  in  front  of  and 
'  behind  the  shock,  respectively,  after  k 

incremental  steps  of  size  h. 

A  potentially  confusing  point  is  raised  by  the  notation  u[t'(kh)]  used 
in  Eq.  (4-2).  When  -we  use  the  Earnshaw  solution  analytically,  -we  compute 
the  value  of  u  as  a  function  of  t'  for  any  point  in  the  (u,t')  plane. 
Note,  however,  that  the  value  of  u  associated  with  a  particular  wavelet 
is  fixed.   In  our  algorithm  we  distort  the  time  base  using  Eq.(3-6), 
calculating  the  new  time  of  arrival  of  some  particular  wavelet.  The 
values  of  u  in  our  application  of  the  Earnshaw  solution  remain  constant. 
Hence,  we  could  conceivably  use  the  source  values  of  u  in  Eq.  (4-2). 
In  our  implementation  of  Eq.  (4-2),  however,  we  shift  and  renumber 
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the  values  of  u  and  t'  associated  with  a  particular  -wavelet.  The  reason 

for  the  shifting  and  renumbering  is  made  clear  below.   The  important 

consideration  here  is  that,  although  the  value  of  u  associated  with 

a  particular  wavelet  is  unchanged  from  its  source  value  by  use  of  the 

Earnshaw  solution,  the  location  of  this  point  in  the  sequence  of  wavelets 

may  be  changed.  Hence,  we  specifically  keep  track  of  individual  values 

of  u  as  a  function  of  kh. 

In  Eq.  (4-2)  note  that  the  value  of  t'[(k+l)h]  is  completely 

determined  by  the  value  of  t'  after  kh  steps  and  the  corresponding  value 

of  u  after  kh  steps.   Similarly,  the  value  of  t  '  after  (k+l)h  steps  is 

dependent  only  upon  the  value  of  t  '  after  kh  steps  and  the  mean  particle 

velocity  at  the  shock  after  kh  steps.   Use  of  Eq.  (4-3)  requires  that  the 

value  of  t  '(kh)  be  retained  for  use  in  the  calculation  of  t  '[(k+l)h]. 
s  s 

In  a  complicated  wave  having  many  shock  fronts,  this  requirement  could 

become  difficult  as  one  shock  merged  with  another.   The  potential 

difficulty  may  be  avoided  by  substituting  for  up  and  u  in  Eq.  (4-3) 

using  Eq.  (4-2)  evaluated  at  t  '(kh).   Carrying  out  the  above  steps, 

s 

we  obtain 

ts'C(k+l)h]  =  tg'(kh)  -  |  jt2'(kh)  -  t2*[(k+l)h]  +  t1'(kh) 

(4-4) 

-  t1'[(k+l)h]l 
Recognizing  that  t  '(kh)  =  t  '(kh)  =  t  '(kh),  we  obtain  the  important 

J_  cL  S 

result, 

t  '[(k+l)h]  +  t.  '[(k+l)h] 
y[(k+l)h]  =-£ g-i .         (4-5) 
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Equation  (  4-5)*  is  a  solution  of  Eq.  (4-3)  and  provides  a  very  simple 
method  for  placing  the  shock.  We  stress  that  Eq.  (4-5)  is  a  solution 
of  Eq.  (4-3)  only  for  a  small  step  size,  that  is,  it  is  a  solution  only 
of  the  finite  difference  form  of  the  shock  occurrence  equation. 

3-   The  Weak-Shock  Theory  Wave  Propagation  Algorithm 

A  computer  algorithm  implementing  Eqs.  (4-2)  and  (4-5)  has 
been  written.   The  computer  language  used  in  this  algorithm  is 
FORTRAN  IV.  We  shall  explain  in  detail  here  the  actual  steps  used  in 
the  algorithm.   Because  the  algorithm  itself  is  rather  complicated,  the 
reader  is  encouraged  to  refer  frequently,  -while  reading  this  section, 
to  the  program  listing  and  flowchart  that  appear  in  Appendix  A. 
Throughout  the  ensuing  discussion  reference  to  specific  variables, 
arrays,  members  of  arrays,  and  so  on  is  in  FORTRAN  notation.   Thus,  the 
Ith  element  of  the  T  array  is  T(l). 

The  algorithm  is  implemented  in  a  subroutine  WAVEPROP  that  is 
called  by  a  driver  program.  We  shall  limit  our  interest  in  the  present 
section  to  this  subroutine.  A  simplified  flowchart  of  WAVEPROP  appears 
in  Fig.  ( 4-l) . 

Before  we  commence  the  detailed  discussion  of  WAVEPROP,  we 
shall  define  and  discuss  a  minor  problem  that  exists  in  WAVEPROP.  We 
call  this  problem  the  "end  point"  problem.   Equation  (4-2)  distorts 
values  in  the  T  array.   Hence,  it  is  possible  for  the  value  of  time  of  a 


Equation  (4-5)  is  a  finite  difference  approximation  of  the  "Landau  Rule" 
The  "Landau  Rule"  is  fully  discussed  by  Landau  and  Lifshitz  (1959). 
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point  close  to  either  end  of  the  data  to  be  distorted  out  of  the  original 
■waveform  time  base.  We  refer  to  a  distortion  of  this  nature  as  an 
"end  point"  problem.  The  reason  it  is  a  problem  is  that  in  writing 
WAVEPROP  we  must  make  some  assumption  concerning  particle  velocity 
values  outside  the  original  time  base.   In  the  present  version  of 
WAVEPROP  we  assume  that  particle  velocity  values  outside  the  original 
time  base  are  unknown.   Hence,  if  the  time  value  of  a  waveform  point  is 
distorted  outside  the  original  time  base  we  are  unable  to  compute  the 
new  waveform  at  that  point.  We  recognize  that  other  assumptions 
concerning  conditions  outside  the  original  time  base  could  have  been 
made  and  implemented.   The  present  version  of  WAVEPROP  could  be  modified 
to  accommodate  these  new  assumptions. 

Subroutine  WAVEPROP  contains  two  self-checking  tests  that 
monitor  for  end  point  problems  and  upon  detection  halt  program  execution. 
The  tests  are  identified  in  WAVEPROP  as  the  variables  TLEFT  and  TRIGHT. 
When  the  program  is  halted  a  diagnostic  message  is  printed,  the  form  of 
which  is* 


EXCEEDED  TIME  BASE  TO 


LEFT 


RIGHT 


DISTANCE  =  YY   KK   ITERATIONS 


* 


We  omit  normal  punctuation  in  lines  containing  FORTRAN  statements 
where  inclusion  of  the  punctuation  could  give  the  impression  that  the 
punctuation  is  part  of  the  FORTRAN  statement. 
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The  word  LEFT  or  RIGHT  indicates  that  some  T  array  value  has  been 
distorted  to  a  value  less  than  the  original  value  of  T(l)*  or  greater 
than  the  original  value  of  T  (NPOINTS),  respectively.  The  value  of 
distance  printed  (YY)  is  the  value  of  the  variable  DISTANCE  in  the  call 
to  WAVEPROP.   The  number  of  iterations  printed  (KK)  is  equal  to  the 
value  of  the  index  being  used  in  loop  1  [see  Fig.  (^-l)]  at  the  point 
of  program  interruption. 

End  point  problems  may  be  easily  avoided.  We  have  used  two 
methods,  depending  upon  the  type  of  wave  that  we  are  considering.   For 
periodic  deterministic  waves ,  we  select  the  phase  of  the  initial  waveform 
in  such  a  fashion  that  when  shocks  form,  they  form  in  the  center  of  the 
cycle,  and  the  particle  velocities  associated  with  T(l)  and  T(NPOINTS) 
are  zero.  When  studying  nonperiodic  deterministic  waveforms  or  noise 
waveforms,  we  place  a  buffer  of  points  having  zero  particle  velocity 
at  both  ends  of  the  original  waveform.   The  subject  waveform  can  then 
expand  its  time  base  during  computer  propagation  into  the  buffer  zones. 
For  the  random  noise  waveforms  that  we  investigate  later  in  this  work, 
buffer  zones  approximately  five  percent  as  long  as  the  total  original 
waveform  were  successfully  used. 

The  subroutine  call  is 

CALL  WAVEPROP(U,  T,  DISTANCE,  DISINC,  NPOINTS) 


In  the  form  of  WAVEPROP  presented  here,  we  are  required  to  set  T(l)=0.0 
However,  T(l)=0.0  is  merely  a  convenience  and  any  value  of  T(l)  is 
acceptable.   If  a  value  of  T(l)  not  equal  to  0.0  is  used,  the  TLEFT 
test  must  be  changed  to  use  this  new  value. 
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Call  WAVEPROP  (U,  T,  DISTANCE,  DISINC,  NPOINTS ] 


Dimension  arrays,  set  constants  and  calculate  the  number 
of  incremental  steps  DISTORS  of  length  DISINC  needed  to 
propagate  the  wave  out  to  the  distance  DISTANCE. 


J  =  1,  NDISTORS 


Set  counting  and  shifting  variables  and  arrays 


■W  2)  1  =  1,  NPOINTS 


The  computations  in  this  loop  sequentially  distorts  each 
point  in  the  T  array  in  accordance  with  equation  (4-2) 
and  detects  shocks  by  comparing  consecutive  elements  of 
the  T  array  for  multivaluedness.  The  times  of  the  shocks 
are  determined  with  the  application  of  equation  (4-5) 


Perform  self-checking  test. 
If  error  is  detected  a  diagnostic  message  is  printed. 


I 


Set  counting  and  shifting  variables  and  arrays 


K  *  1,  NPOINTS 


The  computations  in  this  loop  adjust  particle  velocities 
at  the  shocks  using  a  point-slope  relation.  The  T  array 
is  monitored  for  merging  of  shocks.  When  a  set  of  merging 
shocks  is  detected  it  is  replaced  by  a  single  shock  of 
corrected  amplitude. 


(Return) 


X 

(end) 
FIGURE  4  -  1 


A  Simplified  Flowchart  of  Subroutine  WAVEPROP 
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where 

U        is  a  real  one -dimensional  array  describing  the  particle 
velocity  of  a  point  on  the  subject  -waveform  in  either 
English  or  metric  units.  If  metric  units  are  used,  the 
small-signal  sound  speed  c  must  be  supplied  in  metric 
units.  If  English  units  are  used,  c  at  72°F  is  already 
supplied  in  the  program. 
T        is  a  real  one -dimensional  array  describing  the  time 

coordinate  (in  seconds)  for  a  particular  point  on  the 
subject  waveform.  The  doublet  (T(l),  U(l))  then 
describes  the  Ith  waveform  point. 
DISTANCE  is  a  real  variable.  It  is  the  distance  over  which  we 

wish  to  propagate  the  wave  described  by  U  and  T. 
DISINC    is  a  real  variable.   It  is  the  incremental  step  size 

used  to  cover  the  distance  DISTANCE. 
NPOINTS   is  an  integer  variable.   It  is  the  number  of  waveform 
points . 
The  wave  to  be  propagated  enters  the  subroutine  and  is  returned  through 
the  formal  parameters  U  and  T.   The  quantity  NPOINTS  must  be  entered  as 
a  variable,  not  a  constant,  since  WAVEPROP  changes  NPOINTS  and  returns 
a  new  value .  WAVEPROP  enables  us  to  calculate  the  distortion  of  the 
wave  over  the  distance  given  by  the  real  variable  DISTANCE.   The  number 
of  propagation  steps  NDISTORS  is  set  by  the  relation 

NDISTORS  =  DISTANCE/DISINC  +  0.5  (^-6) 

Computer  operation  causes  NDISTORS  to  round  off  to  the  nearest  integer. 
If  DISTANCE  is  not  an  integral  multiple  of  DISINC,  the  number  of  steps 
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computed  may  "be  in  error  by  one.   The  possibility  of  such  an  error  is 
avoided  by  requiring  that  DISTANCE  be  an  integral  multiple  of  DISINC . 

Referring  to  Fig.  (^-l)  we  observe  that  WAVEPROP  consists  of  an 
outer  DO  loop  (loop  l)  that  is  stepped  through  the  NDISTORS  propagation 
steps  and  two  inner  DO  loops  (loops  2  and  3)  that  are  stepped  through  the 
NPOINTS  points.   For  each  propagation  step  the  entire  waveform  must  be 
distorted  as  prescribed  by  Eq.  (4-2),  the  shocks  (if  any  are  present) 
placed  in  accordance  with  Eq.  (4- 5)j  and  the  velocities  at  the  shocks 
adjusted  to  reflect  the  shock  growth  or  decay.   In  addition,  occurrences 
of  shock  merging  must  be  detected  and  a  single  shock  described  at  the 
merger  location.   In  the  ensuing  discussion  of  WAVEPROP  we  have,  for 
clarity,  described  some  operations  as  occurring  serially.   In  reality 
several  of  these  operations  are  performed  in  a  parallel  fashion. 

The  calculations  required  to  perform  the  tasks  outlined  in 
the  preceding  paragraph  are  divided  between  loops  2  and  J  as  shown  in 
Fig.  (h-l) .      The  function  of  loop  1  is  merely  to  step  the  subject  wave 
through  the  NDISTORS  steps  thereby  propagating  the  wave  the  distance 
DISTANCE.   For  each  NDISTORS  step,  of  course,  loops  2  and  3  are  used 
to  perform  the  necessary  computations. 

Prior  to  each  entry  into  loop  2,  arrays  NS  and  TS  and  variable 
NSHIFT  are  set  equal  to  zero.  The  quantities  TS,  NS  and  NSHIFT  are 
utilized  for  counting  and  shifting  purposes  that  will  become  clear  as 
we  proceed. 

The  first  task  performed  using  loop  2  is  the  sequential 
distortion  of  each  point  in  the  T  array  and  the  detection  of  any 
mult i value dness  in  T  following  the  distortion  calculation.   The 
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distortion,  computed  using  Eq.  (4-2),  is  that  suffered  by  a  particular 
wavelet  in  propagating  a  distance  DISINC .  Mult i value dness  in  T  is 
detected  by  comparing  each  T(l)  to  T(l-l).   If  for  any  value  of  I 
we  find  that 

T(I)  *  T(I-1)    ,  (fc-7) 

multivaluedness  is  indicated.  The  variable  TLOW  is  used  in  the  check 
for  multivaluedness.   Prior  to  entry  into  loop  2,  the  quantity  TLOW 
is  set  as  follows: 

TLOW  =  T(l)  -1.0     .  (4-8) 

The  purpose  of  Eq.  (4-8)  is  to  initially  set  TLOW  out  of  the  time  span 
defined  in  the  T  array.   The  quantity  T(l)  is  then  distorted  and 
compared  to  TLOW.   Since  the  typical  time  shift  per  propagation  step 
is  much  less  than  one  second,  the  initial  setting  of  TLOW  should  preclude 
the  result 

T(l)  <  TLOW    .  (4-9) 

Let  us  assume  that  result  Eq.  (4-9)  is  not  realized.   TLOW  is  then 
changed  to  T(l)  and  compared  to  T(2).   The  scheme  is  repeated,  that  is, 
TLOW  is  set  equal  to  T(l-l)  and  compared  to  T(l)  until  we  obtain  the 
result  (4-7)  or  loop  2  is  completed.   Result  (4-7)  indicates  that 
the  sequence  of  time  values  stored  in  the  T  array  is  no  longer  mono- 
tonic,  which  in  turn  indicates  the  presence  of  a  shock  within  the 
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multivalued  region.*  If  result  (^-7)  is  not  obtained,  no  shocks  are 
present. 

The  location  of  the  shock,  and  the  number  of  points  -within 
the  multivalued  region  must  now  be  determined.  We  make  the  determina- 
tion by  proceeding  as  before,  comparing  TLOW  to  each  T(l).  When  -we 
obtain  the  result 

T(l)  >  TLOW    ,  (1^-10) 

the  end  of  the  multivalued  region  has  been  identified.  In  accordance 
■with  Eq.  (^-5),  the  shock  is  then  located  midway  between  the  vertex** 
points  of  the  multivalued  region.   If  it  is  the  first  shock  found,  it 
is  labeled  NS(l)  and  its  location  is  stored  in  TS(l).   Shocks  detected 
later  are  labeled  NS(2),  NS(3),  etc.,  and  their  locations  stored  in  TS(2), 
TS(3)j  etc.   If  there  are  -waveform  points  having  indices  bet-ween  the 
indices  of  the  t"wo  vertex  points  they  are  counted  and  their  number  stored 
in  NSHIFT.  The  in-between  points  are  now  unnecessary  since  the  region 
is  completely  described  for  our  purposes  by  its  vertex  points.  Hence, 
in  order  to  increase  computational  efficiency  the  NSHIFT  unnecessary 
points  are  shifted  out  of  the  U  and  T  arrays  and  points  having  indices 
greater  than  the  index  of  the  right  vertex  are  renumbered. 


*  We  use  the  term  multivalued  region(s)  to  indicate  those  portions  of 
the  sequence  of  times  stored  in  the  T  array  that  are  not  monotonic. 

**  The  right  and  left  vertex  points  of  a  particular  multivalued  region 
define  the  time  limits  of  the  region,  the  left  vertex  point  being  the 
low  time  limit  and  the  right  vertex  point  the  high  time  limit.   In 
Fig.  (U-2b)  the  point  (i-l)  is  the  right  vertex  and  the  point  I  is  the 
left  vertex  of  the  multivalued  region  pictured. 
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An  illustration  of  the  occurrence  of  multivaluedness  is 
shown  in  Fig.  (k-2) . 


FIGURE  4-2 

A  SECTION  OF  A  FINITE-AMPLITUDE  WAVE  ILLUSTRATING 
THE  COMPUTATIONS  PERFORMED  IN  LOOP  2.  THE  CIRCLES 
ARE  WAVEFORM  POINTS. 

a.  THE  WAVEFORM  AFTER  DISTORTION  BY  EQ.  (4-2). 

b.  THE  WAVEFORM  FOLLOWING  THE  RENUMBERING  OF  THE 
POINTS  AND  THE  PLACEMENT  OF  THE  SHOCK  (TS(1)). 


Figure  (4-2a)  depicts  a  section  of  an  arbitrary  wave  between  the 
points  (l-2)  and  (l+2).   The  multivalued  waveform  in  Fig.  (U-2a)  has 
been  generated  using  Eq.  {k-2).      Proceeding  as  described  above  identifies 
the  (i-l)  point  and  the  (i+l)  point  as  the  vertex  points  of  the  multi- 
valued region.  A  shock  labeled  TS(l) (we  assume  this  is  the  first  such 
region  detected)  has  been  placed  midway  between  the  vertex  points  in 
Fig.  (4-2b),  the  Ith  point  shown  in  Fig.  (^-2a)  has  been  shifted  out, 
and  appropriate  points  have  been  renumbered.   In  Fig.  (4-2b)  the  vertex 
points  are  the  Ith  and  (i-l)th  points.   Figure  (lj-2b)  illustrates  the 
extent  of  the  loop  2  computations.   The  wave  is  still  multivalued  and, 
while  the  shock  has  been  located,  there  has  been  no  adjustment  of  the 
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multivalued  -wave  at  the  shock  location.  These  operations  will  be 
performed  in  loop  3>  where  the  final  -waveform  is  computed. 

After  the  completion  of  loop  2,  we  test  via  program  steps 
in  WAVEPROP,  for  end  point  problems.  We  then  proceed  to  loop  3  unless 
an  end  point  problem  has  been  detected,  in  which  case,  the  program 
is  halted. 

The  propagation  of  the  wave  over  the  distance  DISINC  is 
completed  by  use  of  loop  3-   The  calculations  performed  in  loop  3 
complete  the  waveform  corrections  at  the  shocks  begun  in  loop  2,  detect  # 
shock  mergers,  and  describe  the  waveform  at  the  point  of  shock  merger. 
Recall  that  in  loop  2  there  was  no  adjustment  of  the  magnitude  of  the 
particle  velocity  values.  All  that  we  did  was  to  implement  the  Earnshaw 
solution  and  calculate  the  time  of  occurrence  of  each  shock.   The 
Earnshaw  solution  does  not  affect  the  magnitude  of  the  particle  velocity 

of  the  waveform  points.  Hence,  all  actual  adjustment  of  particle 

velocity  magnitude  at  the  various  waveform  points  is  performed  in 

loop  3» 

The  arrays  TS  and  NS  and  the  variable  NSHIFT  are  again 

zeroed  before  we  start  the  loop  3  computations.* 


The  variable  NSHIFT  is  zeroed  because  the  unnecessary  points  defined 
in  loop  2  have  been  shifted  out  and  zeroing  NSHIFT  is  a  preparatory 
step  for  counting  any  additional  unnecessary  points  defined  during  the 
loop  3  computations.   The  rationale  behind  the  zeroing  of  TS  and  NS  is 
not  quite  so  straightforward.   The  arrays  TS  and  NS  are  used  in  loops  2 
and  3  to  store  shock  locations  and  identifying  numbers.  We  use  the 
information  stored  in  TS  and  NS  within  both  loops  2  and  3-   It  turns 
out  to  be  more  efficient,  however,  to  zero  TS  and  NS  following  the 
loop  2  computations  and  recompute  the  shock  locations  and  identifying 
numbers  during  the  loop  3  computations . 
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Let  us  assume,  as  we  start  the  loop  3  computations,  that  the 
time  sequence  stored  in  the  T  array  contains  one  or  more  multivalued 
regions.   If  it  does  not,  the  results  of  the  loop  3  computations  merely 
confirm  the  monotonicity  of  the  time  sequence  and  control  passes  back 
to  loop  1.   Thus,  assuming  the  presence  of  multivaluedness  in  the  time 
sequence,  we  again  scan  the  T  array  in  order  to  locate  the  various 
multivalued  regions  and  their  vertex  points.   The  scan  of  the  T  array 
is  done  in  exactly  the  same  manner  as  was  done  in  loop  2.   Let  us 
assume  that  we  have  located  the  right  vertex  of  a  particular  multivalued 
region.  We  now  commence  a  backward  scan  through  the  T  array,  searching 
for  the  first  waveform  point  having  a  time  value  less  than  the  time  of 
the  shock  associated  with  the  region.  The  purpose  of  this  backward 
scan  is  to  assist  in  adjusting  the  particle  velocity  in  front  of  the 
shock.   Let  us  call  the  T  array  member  of  the  desired  point  T(j). 
The  points  T(j)  and  T(j+l)  straddle  the  shock  time.   The  value  of 
particle  velocity  immediately  in  front  of  the  shock  is  computed  by  using 
T(j),  T(J+l),  U(j),  U(J+l),  and  the  shock  time  in  a  simple  interpolation 
relation.   In  the  actual  computation  in  loop  3  "the  shock  time  is  called 
TAVG.   The  value  of  particle  velocity  in  front  of  the  shock  is  given  by 
the  relation 

U(J+1)  =  U(J)+  ll^]]    I  f^]    [TAVG-T(J)]    .       (H-U) 

The  T  array  is  then  adjusted  so  that 

T(J+1)  =  TAVG    .  (*+-12) 
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The  net  effect  of  the  use  of  Eqs.  (^-ll)  and  (^-12)  is  to  locate  a 
■waveform  point  immediately  in  front  of  the  shock. 

The  next  computation  performed  after  placement  of  the  "waveform 
point  immediately  in  front  of  the  shock  is  the  removal  of  unnecessary 
points.   Let  us  assume  that  in  the  aforementioned  backward  scan  one  or 
more  -waveform  points  were  tested  before  the  desired  point  in  front  of 
the  shock  was  found.   These  intervening  points  are  now  unnecessary 
because  only  two  points  are  required  to  describe  a  shock.  As  was  done 
in  loop  2,  the  intervening  points  are  shifted  out  of  the  U  and  T  arrays . 

We  now  commence  a  forward  scan  of  the  T  array  starting  at 
the  shock  time.   The  goal  of  forward  scan  is  to  locate  the  first  point 
having  a  time  greater  than  the  shock  time.   Once  we  have  located  the 
desired  point,  the  particle  velocity  and  time  of  the  point  just  behind 
the  shock  are  calculated  using  relations  similar  to  Eqs.  (U-ll)  and  (1+-12). 
Then  any  additional  unnecessary  waveform  points  are  shifted  out. 

Shock  merging,  mentioned  earlier,  is  handled  as  an  integral 
part  of  the  loop  3  calculations  previously  described.  Recall  that  after 
we  initially  locate  each  shock,  the  T  array  is  scanned  backwards  in  order 
to  discover  the  first  waveform  point  in  front  of  the  shock.   During  the 
backward  scan  other  shocks  may  be  encountered  before  the  desired  point 
is  found.   In  that  case,  the  shock  being  considered  will  have  overtaken 
and  passed  the  other  shocks  that  were  encountered.   The  fact  that  other 
shocks  have  been  overtaken  does  not  alter  at  all  the  loop  3  computa- 
tions previously  described.   The  only  additional  computation  required 
is  the  adjustment  of  the  TS  and  NS  arrays  to  reflect  the  adjustment  in 
the  location  and  number  of  shock. 


78 


A  particular  example  of  shock  merging  appears  in 
Fig.  (J+-3)-   In  considering  the  -waveforms  in  Fig.  (^-3),  the  reader 
is  cautioned  that  these  drawings  exaggerate  the  distortion  that  will 
occur  in  one  propagation  step.   In  using  the  difference  Eqs .  (^-2)  and 
(1+-5)  we  assume  that  the  distance  step  size  is  very  small.   Typically, 
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Figure  4-3 
A  section  of  a  finite -amplitude  wave  illustrating  the  computations 
performed  in  loop  3.   (a)    The  wave  upon  entry  into  loop  3.   (b)    The 
wave  has  been  partially  corrected  by  the  formation  of  the  shock 
between  points  2  and  3  and  the  calculation  of  TS(2).   (c)    The  wave 
has  been  further  corrected  by  shifting  out  data  points  in  the  overlap 
region,   (d)    The  final,  fully  corrected,  waveform. 

the  change  in  a  time  coordinate  during  one  propagation  step  of  size  0.1  ft 

is  less  than  10   sec.   In  Fig.  (^+-3)  the  waveform  points  have  "been 

numbered  rather  than  given  the  more  general  labels  used  elsewhere  in 

this  chapter.  Figure  (^-3a)  depicts  the  subject  waveform  upon  entry 

into  loop  3«  Figures  (i+-3b  and  3c)  show  the  waveform  as  we  are  processing 

it  during  the  course  of  loop  3«   In  Fig.  (14-3b)  the  first  shock  TS(l)  has 

been  located,  points  2  and  3  have  been  corrected,  and  the  second  shock 

TS(2)  has  been  located  to  the  left  of  TS(l).   The  fact  that  TS(2)  is  to 

the  left  of  TS(l)  indicates  a  shock  overtaking  situation.  The  first  step 
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in  the  correction  of  this  situation  is  indicated  in  Fig.  (4-3c)  where 
waveform  points  "between  points  2  and  5  in  Fig.  (4-3b)  have  been  shifted 
out  and  points  5  and-  6  renumbered.   The  final  corrected  waveform  is 
shown  in  Fig.  (4-3d);  note  that  the  shock  time  has  been  labeled  TS(l). 
In  the  following  sequence  of  figures,  the  overall  effects  of 
nonlinear  propagation,  as  computed  using  WAVEPROP,  are  illustrated.   For 
simplicity  a  source  signal  consisting  of  straight  line  sections  has  been 
used.   In  Fig.  (4-4)  the  source  waveform  is  shown.   There  are  no  shocks 
present  initially  and  the  signal  contains  a  number  of  local  extrema.   In 
Fig.  (4-5)  the  signal  has  been  propagated  12  ft  using  WAVEPROP;  several 
shocks  have  formed  and  some  of  the  small  shocks  have  already  been  overtaken. 
Figure  (4-6)  shows  the  waveform  at  24  ft.   The  waveform  is  now  approaching 
the  asymptotic  form  described  in  Chapter  3.   The  reader  should  note  that  the 
vertical  scales  in  the  three  figures  are  different. 
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FIGURE  4-4 

THE  SOURCE  WAVEFORM 

WAVEFORM  EXTREMA  HAVE  BEEN  LABELED  TO  ASSIST  THE  READER  IN 
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FIGURE  4-5 
THE  WAVE  SHOWN  IN  Fig.  4-4  PROPAGATED  12  ft 
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FIGURE  4-6 

THE  WAVE  SHOWN  IN  Fig.  4-4  PROPAGATED  24  ft 
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THE  SHOCK  DESCRIBED  BY  THE  SEGMENT  CH 
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*+.   The  Proper  Form  for  Applying  the  Attenuation  and  Dispersion 
Corrections 

In  Chapter  3  we  theorized  that  attenuation  and  dispersion  can 
be  incorporated  into  the  propagation  algorithm  by  assuming  that  small- 
signal  tube  wall  theory  may  be  used  even  for  finite -amplitude  waves 
provided  the  distance  increment  of  propagation  is  small. 

In  this  section  some  mathematical  details  associated  with 
this  incorporation  are  considered.   Nonlinear  distortion,  as  implemented 
using  WAVEPROP,  is  computed  in  the  time  domain.  The  corrections  for 
attenuation  and  dispersion  are  made  in  the  frequency  domain.   Trans- 
formations between  the  time  and  frequency  domain  are  performed  by  use 
of  a  fast  Fourier  transform  (FFT)  computer  subroutine.   The  FFT  is  a 
harmonic  analysis.   Hence,  we  must  modify  Eq.  (3-65)  as  follows: 

-*/na  (l+j)x 

C  '  =  C  e    °         ,  (^-15) 

n     n  ' 

where  a  is  the  value  of  a.   at  the  fundamental  frequency.   In  the  study 
of  random  waveforms  the  fundamental  frequency  is  defined  as  the  reciprocal 
of  the  time  duration  of  the  sample  waveform.   The  actual  method  of 
application  of  Eq.  (4-13)  is  described  below. 

At  any  distance  x  the  waveform  being  propagated  F(t')  can  be 
approximated  by  the  finite  Fourier  series, 

N 

F  ( t ' )  =    y    U     cos  nco  t '    +  UI     sin  nco  t '  ,  ( k-lk) 

Nv      '        ^      n  o  n  o  > 
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■where 

to    =  the  fundamental  frequency, 

F-,(t')  =  N  term  approximation  of  F(t'),    and 

U     and  UI     =  Fourier  coefficients  of  the  nth  series  term, 
n  n 

We  may  express   the   complex  amplitude   of  the  nth  term  of  Eq.    (U-l^)   as 

follows :* 


on  =  un-juin        .  (W5) 


Substituting  Eq.    (i+-15)   into  Eq.    (k-lj) ,   we  obtain 


C    '    =  U    '    -  JUI 
n  n  n 


■where 


-4nOC  x 

U  '  =  e         ° 

n 


U      cos  Vn  a,  x  -   UI     sin  Vn  a  x         , 
L  n  o  n  o  J 


and 


Vna  x 


UI   '   =  e 
n 


°     UI      cos  vn  a  x  +  U      sin  vn  a  x 
L    n  on  o  J 


(^-16) 


The  FFT  used  in  the  actual  computations   calculates  the   complex  conjugate 
of  the  transform.      Hence,   we   shall  use  the  complex  conjugate   in  our 
frequency  domain  analysis. 
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Hence,  the  proper  form  of  F  (f)  after  the  attenuation  and  dispersion 


N 
corrections  is 


C, 


N 


n=l 


^ 


na  x 


U  e       cos  (noo  t'  -  vn  Oi   x) 
n  o        o 


(^17) 


-Jnoc  x 

+  UI  e      sin  (noo  t1  -  vn  a  x) 
n  o        o 


where 


n 

F  (t')  =  the  corrected  form  of  F  (tf ). 

The  value  of  x  used  in  Eqs.  (h-l$),    (h-l6) ,    and  (U-17)  is  not  a  distance 
from  the  origin;  rather  it  is  the  distance  since  the  last  waveform 
correction  for  attenuation  and  dispersion.  In  the  computer  programs 
that  implement  modified  weak-shock  theory  (see  Section  5  of  this  chapter), 
the  selection  of  the  proper  propagation  distance  between  attenuation 
and  dispersion  corrections  is  a  key  factor.  We  have  theorized  that  we 
can  use  small-signal  Kirchhoff  theory  when  it  is  applied  over  a  small 
distance  increment.   In  our  programs  this  small  distance  increment  is 
designated  XL.   In  the  developmental  testing  of  our  algorithm,  various 
values  of  XL  were  tried.   In  general,  we  found  that  XL  must  he  small 
hut  not  so  small  that  we  encounter  the  "resampling"  problem  explained 
in  Section  6  of  this  chapter.  While  we  cannot  offer  any  formulae  or 
rules  for  choosing  the  optimum  value  of  XL,  we  do  give  some  guidelines 
in  Section  6  of  this  chapter. 
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In  the  computer  implementation  of  weak-shock  theory,  the 

corrected  coefficients  U  '  and  UI  '  are  computed  exactly  as  shown  in 

n       n 

Eq.  (4-l6).  We  require  that  the  corrected  -waveform  in  the  time  domain 
be  a  real  function  of  time.   Hence ,  the  correction  terms  are  applied  in 
the  frequency  domain  in  such  a  fashion  that  sequence  (U  ',U  ',  . ..,  U  ') 
is  even  and  the  sequence  (UI  ' ,  UI0',  .  ..,  UI  ']  is  odd  about  the 
(n/2+l)th  point. 

5.   The  Modified  Weak-Shock  Algorithm 

In  this  section  we  shall  show  how  the  various  elements  of 
modified  weak-shock  theory  are  fit  together.   Several  subroutines  are 
used  in  our  algorithm.   Subroutine  WAVEPROP  has  already  been  described. 
The  other  subroutines  are  described  in  this  section, 
a.   The  "Driver"  Programs 

In  a  general  sense  the  implementation  of  modified  weak-shock 
theory  is  exactly  as  indicated  in  Fig.  (3-8) •   The  waveform  to  be 
studied  is  read  into  the  computer  and  propagated  a  short  distance  using 
WAVEPROP.   The  resulting  distorted  waveform  is  Fourier  transformed 
and  the  attenuation  and  dispersion  corrections  are  made.   The  corrected 
frequency  domain  result  is  then  inverse  transformed  and  another  propa- 
gation step  is  taken.   In  this  way,  the  wave  is  propagated  out  to  the 
desired  distance.   In  Appendix  B  typical  driver  program  listings  and 
flowcharts  for  two  cases  of  interest  are  presented.   The  first  of  the 
programs,  PROGRAM  ATENDIS5,  was  used  to  generate  the  theoretical  predic- 
tions for  initially  sinusoidal  waveforms.   The  second,  PROGRAM  NSATDIS, 
was  used  for  noise  waveforms.   The  fundamental  method  in  both  cases  is 
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identical.   The  programs  differ  in  some  minor  peripheral  aspects.  Each 
program  contains  provision  for  propagating  the  subject  wave  by  both 
weak-shock  theory  (i.e.,  nonlinear  effects  only)  and  modified  weak-shock 
theory  (i.e.,  nonlinear  effects  plus  tube  wall  attenuation  and  dispersion 
corrections).   Often  the  program  user  does  not  require  that  both  sets 
of  computations  be  done.   In  this  case,  the  undesired  portion  of  the 
program  may  easily  be  bypassed  by  the  insertion  of  appropriate  program 
statements .   The  variables  used  in  the  two  programs  are  defined  in  the 
listings  which  appear  in  Appendix  B.   Typical  values  of  some  of  these 
variables  are  discussed  in  Section  7  of  this  chapter. 

Early  in  the  listings  of  both  NSATDIS  and  ATENDIS5,  the 
user  is  cautioned  that  the  dimension  of  U  and  T  must  exceed  by  two  the 
number  of  points  needed  to  describe  the  wave.   The  reason  for  the  extra 
points  is  simple  although  somewhat  obscure.  NSATDIS  and  ATENDIS  both 
call  a  plotting  routine  called  WAVEPLOT.  WAVEPLOT  in  turn  calls  a 
scaling  routine  (SCALE)  that  scales  numbers  to  values  appropriate  for 
the  plot  size  desired.  We  utilize  the  two  extra  locations  in  the  U  and  T 
arrays  to  store  l)  the  array  minimum  value,  and  2)  the  incremental  point 
separation,  both  of  which  are  calculated  by  use  of  subroutine  SCALE, 
b.   The  Subprograms 

Five  subprograms  called  by  programs  ATENDIS5  and  NSATDIS  are 
described  in  this  subsection.  Unless  otherwise  indicated,  flowcharts 
and  listings  for  these  subprograms  appear  in  Appendix  C.   The  functions 
of  three  of  the  subprograms  are  obvious .   The  functions  of  subprograms 
RESAMPLE  and  FIXER,  however,  are  not  obvious.   In  Section  6  we  include 
a  detailed  discussion  of  the  functions  of  RESAMPLE  and  FIXER. 
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1)  Function  ALPHA.   This  function  is  used  to  calculate  the 
value  of  the  Kirchhoff  attenuation  and  dispersion  coefficient  a.  We 
implement  Eq.  (3-6*0  in  function  ALPHA.   The  function  call  is 

A  =  ALPHA(N,F0) 

■where 

N  =  the  harmonic  number,  and 

FO  =  the  fundamental  frequency. 

A  =  a  storage  location  for  the  number  returned  by  ALPHA. 
Since  ALPHA  is  so  simple,  neither  a  flowchart  nor  listing  of  ALPHA,  has 
been  included.   The  radius,  R,  must  be  provided  by  the  user  in  ALPHA. 

2)  Subroutine  FFTC .   This  subroutine  is  a  fast  Fourier 
transform  (FFT)  available  as  a  library  routine  on  the  CDC  3200  facility 
at  the  Applied  Research  Laboratories,  The  University  of  Texas  at  Austin. 
It  was  written  in  COMPASS  by  James  D.  Beard.   The  subroutine  call  is 

CALL  FFTC(X,Y,M) 

where 

X  is  a  real  one -dimensional  array  containing  the  real  part  of 

the  input  and  output  functions, 
Y  is  a  real  one -dimensional  array  containing  the  imaginary  part  of 

the  input  and  output  functions, 
M  The  number  of  points  in  x  and  y  to  be  transformed  is  Z  . 
Since  subroutine  FFTC  is  a  standard  library  routine,  its  flowchart 
and  listing  are  not  included  here. 
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3)  Subroutine  RESAMPLE.   This  subroutine  is  used  to 
resample  the  subject  -waveform.   The  subroutine  call  is 

CALL  RESAMPLE  (U,  T,  N,  NPOINTS,  DT) 
where 

U   is  a  real  one -dimensional  particle  velocity  input-output  array, 
T   is  a  real  one -dimensional  time  input -output  array, 
N   is  an  integer  variable  equal  to  the  number  of  output  data 
points  desired, 
NPOINTS  is  an  integer  variable  equal  to  the  number  of  input  data  points, 
DT   is  a  real  variable  equal  to  the  sampling  interval  of  the 
output  data  points. 
The  N  output  data  points  are  computed  by  an  interpolation 
scheme.  A  test  in  the  program  halts  execution  if  interpolation 
proceeds  beyond  the  end  of  the  input  data.   In  this  case,  the  diagnostic 
message 

NN  EXCEEDED  NPOINTS 

is  printed.   The  routine  will  -work  correctly  for  N  both  greater  than 
and  less  than  NPOINTS.   Subroutine  RESAMPLE  requires  a  scratch  space 
of  1200  floating  point  words.   If  the  dimension  of  U  and  T  -were  to 
exceed  1200,  this  scratch  space  -would  have  to  be  increased. 
h)      Subroutine  FIXER.   The  subroutine  call  is 

CALL  FIXER  (U,  T,  PHI,  N,  ISTEP,  PER) 

-where 

U   is  a  real  one -dimensional  particle  velocity  input -output  array, 
T   is  a  real  one -dimensional  time  input-output  array, 
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PHI  is  a  real  one-dimensional  output  array,  whose  Ith  element 

contains  the  phase  shift  detected  in  the  Ith  call  to  FIXER, 
N   is  an  integer  variable,  equal  to  the  number  of  points  in  the 
U  and  T  array  plus  1, 
ISTEP  is  an  integer  variable  that  counts  the  calls  to  FIXER, 
PER  is  a  real  variable  equal  to  the  period  in  seconds  of  the 
■waveform  described  by  the  U  and  T  arrays. 

FIXER  does  not  contain  any  diagnostic  tests.   Subroutine  FIXER  calls 
subroutine  RESAMPLE.   Hence,  FIXER  cannot  be  used  -without  RESAMPLE. 
FIXER  requires  2  scratch  arrays  of  dimension  513  floating  point  words 
each.   If  the  dimension  of  U  and  T  "were  to  exceed  515>  these  scratch 
arrays  would  have  to  be  increased. 

5)   Subroutine  WAVEPLOT.   Subroutine  WAVFPLOT  is  used  to 
plot  waveforms.   The  subroutine  call  is 

CALL  WAVEPLOT(U,  T,  DISTANCE,  NPOINTS,  ISW,  TAXIS,  UAXIS,  VA) 

The  definitions  of  the  formal  parameters  are  as  follows : 

U    is  a  real  one -dimensional  particle  velocity  array. 
T    is  a  real  one -dimensional  time  array. 
DISTANCE  is  a  real  variable  equal  to  the  propagation  distance  at 
which  the  plot  is  drawn. 
NPOINTS  is  an  integer  variable  equal  to  the  number  of  data  points 
to  be  plotted. 
ISW    is  an  integer  variable  switch  whose  value  may  be  0  or  1. 
ISW  sets  the  plot  vertical  size.   If  ISW=0  every  call  to 
WAVEPLOT  produces  a  plot  UAXIS  inches  high.   In  subsequent 
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calls  the  vertical  size  of  the  plotted  curve  is  adjusted  to 
maintain  the  proper  amplitude  proportionality  between  the 
waveforms  at  the  first  call  distance  and  at  the  later  distance, 
TAXIS   is  a  real  variable  equal  to  the  length  in  inches  of  the 

horizontal  (time)  axis. 
UAXIS   is  a  real  variable  equal  to  the  length  in  inches  of  the 

vertical  (particle  velocity)  axis. 
VA      is  a  real  variable  that  adjusts  the  vertical  axis  for 

oscilloscope  vertical  amplification.  When  used  with  ISW=1  it 
allows  direct  comparison  between  oscillograms  and  computer 
plots.   If  ISW=0  set  VA=1. 

WAVEPLOT  calls  a  number  of  standard  plotting  subroutines . 
The  plotter  utilized  was  a  CALCOMP  model  565.   The  plotting  subroutines 
are  fully  explained  in  the  Applied  Research  Laboratories  plot  manual 
(1973)-   The  plotting  subroutines  used  were 

a)  PLOT, 

b)  SCALE, 

c)  LINE, 

d)  AXIS, 

e)  PLOTLINE. 


One  additional  plotting  subroutine  called  subroutine  PLOTBUFF 
was  used.  Use  of  PLOTBUFF  allows  us  to  CALCOMP  plot  while  the  computer 
proceeds  with  other  tasks  in  the  program  being  processed.   Parallel 
operation  of  the  plotter  and  computer  saves  considerable  computer  time 
when  several  plots  are  constructed.   Omission  of  PLOTBUFF  slows  program 
execution  but  does  not  otherwise  affect  it. 
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6.   Computational  Details 

In  the  preceding  sections,  the  application  of  the  overall 
modified  weak-shock  algorithm  was  explained.   During  this  explanation 
we  ignored  several  minor  problems  and  their  solutions  contained  in  the 
overall  method.   In  addition  we  have  not,  as  yet,  indicated  how  we 
select  several  important  program  parameters.   These  omissions  are 
rectified  in  this  section, 
a.   Problem  Areas 

1)  "End  Point"  Problem.  This  problem  and  its  solution  have 
already  been  discussed  in  detail  in  Section  3.  We  list  the  end  point 
problem  here  so  that  the  problem  area  subsection  may  be  complete. 

2)  Resampling.   The  sampled  data  used  to  represent  a 
particular  signal  is  equally  spaced  when  it  is  originally  inserted 
into  one  of  the  modified  weak-shock  programs.  A  characteristic  of 
subroutine  WAVEPROP,  however,  is  to  distort  the  time  values  of  sample 
points  by  an  amount  proportional  to  their  particle  velocities. 
There  is  an  additional  adjustment  of  the  waveform  points  at  the  shock 
fronts.   The  end  result  on  the  waveform  points  is  twofold.   The  number 
of  points  is  reduced  as  unnecessary  points  at  the  shock  fronts  are 
shifted  out,  and  the  remaining  data  points  are  no  longer  equally  spaced. 
Unequal  spacing  does  not  cause  any  problems  so  long  as  frequency  domain 
computations  are  not  required.   However,  the  attenuation  and  dispersion 
corrections  are  performed  in  the  frequency  domain,  and  equally 

spaced  waveform  points  are  required  for  use  in  the  subroutine  FFTC. 
Thus  adjustment  of  the  data  points  is  necessary.   In  subroutine 
RESAMPLE  a  point-slope  formula  is  used  to  adjust  the  data  points. 
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A  flow  chart  and  listing  of  this  subroutine  appears  in  Appendix  C . 
Subroutine  RESAMPLE  resamples  the  -waveform  and  returns  the  required 
number  of  equally  spaced  samples.   Since  subroutine  RESAMPLE  is  usually 
followed  by  subroutine  FFTC,  the  number  of  sample  points  produced  by 
RESAMPLE  is  normally  2" +1,  where  n  is  a  positive  nonzero  integer.* 
Application  of  RESAMPLE  cures  the  unequal  spacing  problem  but  at  the 
expense  of  introducing  another  problem:  RESAMPLE  effectively  "unshocks" 
any  shocks  present  in  the  waveform  being  resampled.   The  output  of 
RESAMPLE  consists  of  single  valued,  equally  spaced  waveform  points. 
Hence,  the  minimum  shock  rise  time  after  resampling  is  equal  to  the 
value  of  the  variable  DT  used  in  subroutine  RESAMPLE.  We  call  the 
problem  that  we  have  just  described  the  "resampling"  problem. 
Resampling  is  a  problem  because  the  resampled  waveform  must  be  propa- 
gated using  WAVEPROP  some  short  distance  before  the  shocks  reform  and 
waveform  attenuation  at  the  shock  fronts  resumes .  We  can  minimize  the 
effects  of  resampling  by  minimizing  our  calls  to  RESAMPLE.   However, 
calls  to  RESAMPLE  precede,  and  are  a  necessary  part  of  the  corrections 
for  attenuation  and  dispersion  and  we  have  already  required  frequent 
corrections  (see  Section  h) .     We  have  not  devised  a  method  of  eliminating 
the  resampling  problem.   The  effects  of  resampling  can,  however,  be 
minimized  subject  to  the  Constraint  that  we  make  frequent  Kirchhoff 
theory  corrections  by  careful  selection  of  the  variable  XL.   Our 
guidelines  for  the  selection  of  XL  appear  in  Subsection  6  of  this 
section. 


Subroutine  FFTC  uses  2  sample  points.   The  one  extra  point  produced  by 
RESAMPLE  is  the  first  point  of  the  next  cycle  of  the  wave  and,  since 
it  is  the  first  point,  it  has  zero  particle  velocity. 


3)  Phasing.  The  problem  of  "phasing"  arises  in  the  correction 
of  deterministic  waveforms  for  thermoviscous  attenuation  and  dispersion. 
The  meaning  of  the  term  phasing  is  best  illustrated  by  a  simple  example. 
Consider  a  deterministic  -waveform  that  has  been  distorted  using 
WAVEPROP  and  modified  further  by  the  application  of  correction  factors 
in  the  frequency  domain.  Before  entry  into  the  frequency  domain,  the 
sample  record  consists  of  one  or  more  complete  cycles.   The  record  end 
points  have  zero  particle  velocity.   The  dispersion  correction  in  the 
frequency  domain  changes  the  overall  phase  of  the  reconstituted  time 
•waveform.   The  corrected  -waveform  still  consists  of  one  or  more  complete 
cycles ,  but  the  waveform  end  points  no  longer  have  zero  particle  velocity. 
Hence,  in  the  next  application  of  subroutine  WAVEPROP,  we  will  encounter 
an  end  point  problem. 

Use  of  subroutine  FIXER  enables  us  to  avoid  the  end  point 
problem.   In  FIXER  -we  utilize  the  fact  that  a  full  cycle  of  the  wave 
is  present.   The  end  points  in  the  wave  are  rearranged  to  form  a  full 
cycle  starting  and  stopping  on  a  zero  particle  velocity  point.  The 
phase  difference  between  this  waveform  and  the  input  waveform  is  calcu- 
lated.  The  phase  difference  calculation  is  a  measure  of  the  total 
phase  shift  due  to  dispersion.   For  each  entry  into  FIXER,  the  phase 
shift  is  stored  in  the  appropriate  element  of  the  array  PHI.  A  flow 
chart  and  listing  of  subroutine  FIXER  appears  in  Appendix  C . 

FIXER  is  only  used  when  studying  periodic  deterministic  waves . 
The  same  type  of  phase  shift  difficulty  arises  with  noise  waveforms, 
but  with  noise  signals  another  remedy  is  used.  We  remarked  earlier 
that  noise  signals  are  buffered  at  either  end  by  a  number  of  waveform 
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points  having  zero  particle  velocity.  After  each  dispersion  correction, 
we  merely  rezero  the  first  few  points  in  the  left  buffer  and  last  few 
points  in  the  right  buffer.  The  rezeroing  prevents  the  occurrence  of 
end  point  problems . 

b.   Parameter  Selection 

Programs  NSATDIS  and  ATENDIS5  contain  a  number  of  parameters 
whose  values  are  user  controlled.  Values  for  most  of  the  parameters  are 
not  critical.   In  this  subsection  we  consider  four  key  parameters  whose 
value  significantly  affects  the  results  obtained  using  the  modified 
weak-shock  algorithm. 

1)  SF.   This  parameter  is  a  scale  factor  that  converts  the  input 
data  to  units  of  particle  velocity.  Most  of  our  input  data  was  measured 
directly  from  oscillograms.   The  value  of  SF  contains  the  receiving 
system  sensitivity.   In  programs  NSATDIS  and  ATENDIS5  the  units  of 

SF  are  ft/sec-in. 

2)  XINC.   This  parameter  is  the  incremental  step  size  used  in 
WAVEPROP.   The  value  of  XINC  used  here  was  determined  by  comparing 
results  obtained  using  WAVEPROP  to  results  obtained  using  an  exact 
expression.   Specifically,  we  calculated  the  distortion  of  a  sinusoidal 
wave  as  a  function  of  distance  using  WAVEPROP  and  then  computed  the 
fundamental  through  tenth  harmonic  amplitudes  using  FFTC .  We  compared 
these  results  at  each  distance  at  which  they  were  computed  to  the  funda- 
mental through  tenth  harmonics  calculated  using  Blacks tock's  (l961+) 
exact  equation.  We  found  that  for  values  of  XINC  less  than  one-tenth 

of  the  shock  formation  distance  the  differences  between  the  two  computa- 
tions were  very  small  (less  than  10$  for  all  harmonics  investigated), 
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and  did  not  grow  with  distance.   For  values  of  XINC  greater  than  one- 
tenth  of  the  shock  formation  distance  some  of  the  differences  became 
greater  than  10$.  Using  the  results  of  this  analysis,  we  decided  that 
a  value  of  XINC  equal  to  one -tenth  of  the  shock  formation  distance  was 
appropriate. 

For  the  strongest  sinusoidal  waves  considered  here,  the 
shock  formation  distance  is  3-8  ft  (see  Table  V) .   In  the  study  of 
noise  waves  some  portions  of  the  waves  are  so  strong  that  shocks  form 
very  close  (perhaps  as  close  as  1  ft)  to  the  source  measurement  location. 
Therefore  we  selected  0.1  ft  as  the  value  of  XINC  in  all  analyses 
reported  here.   Cook  (1962)  selected  a  step  size  of  one-twentieth  of 
the  shock  formation  distance.   Our  analysis  indicates  that  one-tenth 
is  adequate. 

When  considering  waves  whose  shock  formation  distance  is 
greater  than  1  ft  larger  values  of  XINC  may  be  used.  We  have  not 
chosen  to  make  this  adjustment  in  the  work  reported  here. 

3)  XL.   This  parameter  is  the  incremental  propagation  distance 
between  successive  attenuation  and  dispersion  corrections.   Its  selection 
is  a  compromise  between  the  necessity  for  frequent  waveform  correction 
and  the  necessity  to  minimize  the  resampling  problem  (see  Subsection  6. a. (2) 
above).  Our  guidelines  for  selecting  XL  are  as  follows: 

1)  Inspect  the  input  record  and  select  the  highest  peak  u  and 

a 

the  lowest  peak  u,  that  may  eventually  become  part  of  a  shock. 

2)  Using  Eq.  (^-2)  calculate  the  net  time  distortion  of  points 
having  particle  velocities  u  and  u,  when  propagated  a 
distance  XL. 
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3)   Calculate  the  point  spacing  of  the  -waveform  at  the  source 
(i.e.,  before  any  distortion). 

k)       Form  ratios  between  the  time  values  computed  in  step  2  and 
the  time  value  computed  in  step  3« 

5)   A  useful  rule  of  thumb  is  that  XL  should  be  selected  in  such 
a  way  that  the  ratios  formed  in  step  h   fall  between  2  and  7- 
The  above  procedure,  while  useful,  is  far  from  perfect  and  some 
experimentation  with  the  value  of  XL  may  be  required.   In  the  results 
reported  in  this  work,  XL  was  selected  to  fall  within  the  limits 
specified  in  step  5. 

k)       M.   The  selection  of  NN,  the  number  of  initial  waveform  points, 
is  a  compromise  between  conflicting  demands.  We  have  already  discussed 
the  problem  associated  with  resampling.   In  order  to  minimize  the 
resampling  problem  while  still  using  small  distance  increments  a  large 
number  of  waveform  points  should  be  used.  A  second  reason  for  selecting 
a  large  number  of  data  points  is  to  avoid  "aliasing".*  Yet  as  the 
number  of  points  increases,  computation  time  also  increases.   The  major 
time  user  in  both  the  ATENDIS5  and  NSATDIS  programs  is  the  subroutine  FFTC 
As  Welch  (1967)  points  out,  the  time  required  to  perform  a  finite  Fourier 
transform  on  a  sequence  of  length  L  is  approximately  KL  log?L,  where  K 
is  a  constant  which  depends  on  the  program  and  the  computer  used. 
Therefore  from  a  cost  viewpoint  we  wanted  as  few  data  points  as  possible. 


"Aliasing"  is  an  effect  that  is  present  in  all  spectral  computation 
schemes,  both  analog  and  digital.   It  is  the  folding  of  frequency  compo- 
nents above  the  Nyquist  frequency  back  into  the  true  spectrum.  Aliasing 
can  be  avoided  if  the  sampling  frequency  exceeds  by  a  factor  of  2 
the  highest  frequency  in  the  sampled  data.  An  excellent  elementary 
discussion  of  the  principles  of  data  sampling,  including  aliasing, 
appears  in  a  book  by  Stein  and  Jones  (1967).   Stein  and  Jones  refer  to 
aliasing  as  a  spectral  ambiguity. 
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Our  solution  to  this  problem  was  obtained  by  trial  and  error.   In  the 
developmental  testing  of  our  algorithm  we  utilized  waves  having 
fundamental  frequencies  between  500  Hz  and  3500  Hz .  We  found  that  by- 
using  512  data  points  per  cycle  we  obtained  excellent  results. 
Increasing  the  value  of  NN  to  1025  points*  gave  negligible  improvement. 
Decreasing  NN  to  257  points  resulted  in  noticeable  loss  of  agreement 
between  our  computed  and  observed  waveforms,   in  the  study  of  noise 
waveforms  the  time  duration  of  the  signals  was  longer  than  in  the  study 
of  sinusoids.   Hence,  more  waveform  points  were  needed.  Again  the  best 
method  of  determining  what  value  of  NN  to  use  proved  to  be  trial  and  error. 
We  usually  started  with  NN=513  points  and  calculated  the  noise  waveform 
after  propagation.  We  then  repeated  the  calculation  at  1025  points  and 
so  on  for  higher  values  of  Z  +1   points .  When  no  appreciable  change 
occurred  between  runs,  it  was  assumed  that  the  number  of  data  points 
had  become  adequate.   For  the  noise  waveforms  reported  in  this  work 
1025  points  appeared  to  be  enough.   The  number  of  data  points  used 
The  number  of  data  points  used  establishes  the  Nyquist 
frequency.   Since  the  reviewing  system  used  in  our  experiments  is  useful 
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*  T\T  N 

In  Program  ATENDIS5  we  use  NN^+i  points.   The  first  2  point: 
describe  a  complete  cycle,  the  last  point  is  the  first  point  of  the 
next  cycle.   The  particle  velocity  associated  with  this  extra  point 
is  zero.   The  function  of  the  extra  point  is  to  prevent  the  occurrence 
of  an  end  point  problem. 

The  frequency  characteristic  of  the  microphones  used  in  our 
experiment  is  constant  up  to  70  kHz.   Because  the  microphones  are 
installed  in  such  a  way  that  the  signal  is  observed  at  a  grazing 
incidence,  we  must  consider  the  directivity  of  the  microphones.  A 
full  discussion  of  the  directivity  pattern  appears  in  Chapter  5.   It 
turns  out  that  the  microphone  directivity  limits  the  bandwidth  of  the 
microphones  in  our  experiment  to  about  27.5  kHz. 
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sampling  rate  commensurate  -with  a  bandwidth  of  27.5  kHz.   Hence  we 
desire  a  55  kHz  sampling  rate.   The  required  number  of  samples  then  is 
55  samples  per  millisecond.   During  the  investigation  of  sinusoids,  we 
considered  a  single  cycle.   The  lowest  frequency  considered  was  500  Hz. 
For  a  500  Hz  wave  110  samples/cycle  are  required  to  meet  the  55  kHz 
sampling  criterion.  We  need  fewer  samples/cycle  for  higher  frequency 
waves .   The  number  of  samples  used  in  subroutine  FFTC  must  be  an  integral 
power  of  two.   The  use  of  512  waveform  points/cycle  easily  meets  all 
requirements  for  the  sinusoidal  waves  investigated  here. 

In  the  study  of  noise  waveforms  the  longest  sample  record 
used  was  16.27  msec.   In  this  case  we  used  NN=1025.   Therefore  the 
sample  rate  was  63  kHz . 

Hence,  we  conclude  that  in  the  cases  investigated  here  the  use 
of  513  sample  points  for  sinusoids  and  1025  sample  points  for  noise 
waves  was  adequate. 
7.   Calculation  of  Power  Spectra 

The  primary  purpose  of  this  whole  research  effort  has  been 
to  investigate  the  change  in  waveform  and  spectral  composition  of 
finite -amplitude  noise  as  it  propagates.  We  decided  that  the  most 
meaningful  spectral  calculations  would  be  those  showing  the  change  in 
the  source  spectrum  due  to  weak-shock  effects  only.   A  computation 
scheme  was  devised  to  calculate  the  power  spectrum  at  various  distances 
from  a  source  of  unbounded  plane  waves.   In  this  program  we  consider 
weak-shock  effects,  but  not  tube  wall  attenuation  and  dispersion.  A 
secondary  reason  for  considering  only  the  weak-shock  effects  is  economic. 
The  application  of  the  tube  wall  effect  corrections  represents  a  major 
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up  to  about  27.5  kHz,-**  -we  -want  to  provide  a  time  portion  of  each  computer 
run.   Hence  in  neglecting  these  corrections  "we  obtain  a  considerable 
saving.   A  schematic  diagram  of  the  program,  program  D0ITALL3,  appears 
in  Fig.  (4-7).  A  detailed  flowchart  and  listing  appears  in  Appendix  D. 
All  arrays  and  variables  used  are  defined  in  the  listing. 

The  heart  of  program  D0ITALL3  is  the  calculation  of  the 
■waveform  distortion  due  to  nonlinear  propagation  and  the  subsequent 
calculation  of  power  spectral  density.  We  minimize  the  effects  of 
resampling  by  not  resampling  the  wave  until  it  has  been  propagated 
the  required  distance.   A  measure  of  economy  in  the  use  of  core  memory 
is  achieved  by  the  multiple  purpose  use  of  arrays.   The  program  has  a 
built  in  capacity  for  ensemble  averaging  power  estimates.  A  scratch 
tape  is  used  to  temporarily  store  computed  spectra  thereby  further 
decreasing  core  requirements.  The  spectral  computation  scheme  is  an 
updated  version  of  the  Blackman  and  Tukey  (1959)  indirect  method. 
The  method  is  indirect  (as  opposed  to  direct)  because  the  autocorrelation 
function  is  calculated  first,  then  the  power  spectral  density  is  computed. 
The  steps  in  the  indirect  method  are  as  follows : 

l)   Let  {X  ,  Xp,  . ..,  X  }  be  a  sequence  of  data  points  and  At  the 
time  interval  between  adjacent  points.   Compute  the  autocorrelation 
function  of  the  sequence  {X  ,  . . . ,  X  } .   The  autocorrelation  function 
is  computed  by  use  of  an  FFT.   The  autocorrelation  function  is  periodic, 
values  computed  using  negative  lags  being  stored  in  the  second  half  of 
the  array. 
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AVERAGE  SPECTRUM  FOR 
EACH  DESIRED  DISTANCE 


FIGURE  4-7 
SCHEMATIC  OF  PROGRAM  DOITALL3 


ARL  -  UT 
AS -73 -596 
FMP  -  DR 
6-26-73 
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2)  Construct  a  sequence  of  mean  lagged  products*  (C  ,  C  ,  . ..,  C  , 
C  _,  .  ..,  C   )  made  up  of  the  first  m  elements  and  the  last  m  elements 

of  the  autocorrelation  function.   The  sequence  C. ,  ....  C„  is  an  even 

1       dm 

function  about  the  (m+l)th  point.  The  time  interval  between  elements  of 
the  sequence  {C  ,  .  ..,  C„  )  is  At,  where  AT=hAt,  h  being  a  positive 
nonzero  integer. 

3)  The  sequence  of  mean  lagged  products  is  then  transformed. 

Since  the  mean  lagged  products  form  the  sample  values  of  an  even  function, 

the  complex  FFT  is  equivalent  to  Blackman  and  Tukey's  cosine  transform. 

The  sequence  of  transformed  values  (Vn ,  V„,  ....  W,  }  is  itself  an 

12'    7   2m 

even  function  symmetric  about  every  integral  multiple  of  m+1.   Blackman 
and  Tukey  (1959)  call  this  sequence  the  sequence  of  raw  spectral  density 
estimates . 

h)       The  refined  spectral  density  estimates  are  then  computed. 
This  computation  is  made  to  combat  the  effects  of  "spectral  leakage" 
Spectral  leakage  is  caused  by  energy  at  frequencies  in  between  the 
computed  harmonics  that  can  be  partially  added  to  the  energy  at  a  computed 
harmonic  because  the  true  spectrum  convolves  with  the  spectrum  of  the 
data  window.   Careful  selection  of  the  data  window  shape  can  greatly 


At  the  time  that  the  book  by  Blackman  and  Tukey  (1959)  "was  written,  the 
mean  lagged  products  were  computed  directly  by  shifting,  multiplying, 
and  adding.   Now,  of  course,  as  Bingham,  Godfrey,  and  Tukey  (1967)  have 
pointed  out,  the  efficient  way  to  perform  this  calculation  is  by  use 
of  the  FFT  and  the  transform  pair  property  of  the  autocorrelation 
function  and  power  spectral  density. 
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reduce  this  leakage.   The  window  used  in  this  work  is  a  hanning 
window, 


D(t)  =|fl+  cos  |21j    |T|  <Tm 


(*-l8) 

=  0  |t|  >  T 


m 


where 


T  =  the  lag  time  of  the  mth  lag  product. 


The  refined  spectral  density  estimates  produced  by  the  indirect 
method  have  equivalent  degrees  of  freedom  in  the  CHI-squared  sense  given 
by  the  relation 


k  = 
r 


nAt  |  1 
mAT   3 


r  =  2,  3,  . . . ,  m 

(^-19) 
half  as  much    r  =  1,  m  ♦  1   , 


where 


k  =  degrees  of  freedom. 
r 

Increased  statistical  stability  is,  of  course,  the  advantage  of  the 

indirect  method.   In  examples  presented  in  Chapter  6  of  this  work 

n=102k,   h=l,  and  m= 128.   Substituting  these  values  into  Eq.  (k-iy) , 
we  obtain 


k  =  7-77    r  =  2,  ...,  m, 
r  (i^-20) 


=  3.88    r  =  1,  m  +  1 
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The  r-1  and  m+1  points  are  the  d.c.  value  and  the  value  at  the  Nyquist 

frequency,  respectively.  The  indirect  method  in  this  case  gives  results 

as  statistically  stable  as  the  ensemble  average  of  3  sets  of  spectral 

estimates  calculated  by  the  direct  method.   In  the  direct  method,  the 

sequence  {X  ,  . . . ,   X  )  is  transformed  using  the  FFT  to  get  the  complex 
1       n 

coefficients,  A.,  .B  .   These  coefficients  are  then  hanned  and  power 

2   2 
terms  of  the  form  A  +B   are  calculated.   The  direct  method,  -while 

computationally  simple,  has  only  2  degrees  of  freedom.   It  would  appear 

that  the  indirect  method,  because  of  its  greater  statistical  stability, 

would  be  preferable.   This  is  not  always  the  case  because  frequency 
resolution  is  lost  when  the  indirect  method  is  used.   For  example,  in 
using  the  sequence  of  mean  lagged  products  we  reduce  the  number  of  data 
points  from  n  to  2m.   Typically,  this  might  be  a  reduction  by  a 
factor  of  k.      The  frequency  resolution  is  degraded  by  this  same  factor. 
Hence,  the  selection  of  the  parameters  used  in  D0ITALL3  dictates  the 
resolution  and  stability  of  the  resulting  power  density  estimates. 
The  theoretical  results  obtained  using  this  method  and  experimentally 
obtained  spectra  are  presented  in  Chapter  6. 


CHAPTER  5  EXPERIMENTAL  EQUIPMENT  AND  PROCEDURES 

1.  Introduction 

In  this  chapter  we  concentrate  on  the  actual  pieces  of 
equipment  and  procedures  used  to  make  the  experimental  measurements 
that  are  reported  in  Chapter  6.   The  second  section  contains  listings 
of  the  various  pieces  of  equipment  used  in  our  experiments.   Within  the 
listings  are  detailed  discussions  of  the  equipment  calibration,  charac- 
teristics, and,  where  applicable,  accuracy  that  directly  affected  our 
experimental  results.   In  the  third  section  we  briefly  outline  each  of 
our  experimental  procedures  and  indicate  the  particular  pieces  of  equip- 
ment from  the  lists  in  Section  2  that  were  used  in  each  procedure.   In 
certain  cases  where  two  or  more  similar  devices  were  available  (e.g., 
two  different  oscillators)  we  indicate  why  the  particular  device  that  we 
used  was  selected. 

The  assumptions  made  in  the  theoretical  development  presented 
in  Chapter  3  place  physical  restrictions  upon  our  experiments.   In 
Section  h   of  this  chapter  these  restrictions  are  analyzed.   The  final 
section  contains  an  analysis  of  the  overall  accuracy  of  our  measure- 
ment equipment. 

2.  Experimental  Apparatus 

A  generalized  schematic  of  the  experimental  arrangement  appears 
in  Fig.  (5-1 ).   All  equipment  was  calibrated,  if  necessary,  before  its 
use.   The  accuracy  of  the  various  calibrations  and  the  precision  with 
which  experimental  data  was  gathered  contribute  directly  to  the  accuracy 
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FIGURE  5-1 

A  GENERALIZED  SCHEMATIC  OF  THE  EXPERIMENTAL 

EQUIPMENT   ARRANGEMENT 


of  the  experimental  results  reported  in  Chapter  6.   Commercially 
available  equipment,  calibrated  as  directed  in  the  applicable  technical 
manual,  in  all  cases  met  or  bettered  the  manufacturer's  nominal  specifi- 
cations.  The  calibration  of  locally  prepared  equipment  will  be  reviewed 
in  this  section.   The  equipment  arrangement  can  be  conveniently  subdivided 
into  a  transmitting  system,  the  pipe,  and  a  receiving  system, 
a.   The  Transmitting  System 

1.  General  Radio  Oscillator  Type  13IO-A. 

2.  Bruel  and  Kjaer  Oscillator  Type  1022. 

3.  Beckman  Universal  Eput  and  Timer  Counter  Model  73^0  AR. 

4.  Krohn-Hite  Filter  Model  315-ACR.   The  function  of  this 
filter  was  to  shape  the  input  spectrum  to  the  sound  source  in  such  a 
fashion  that  energy  outside  the  sound  source  passband  was  not  supplied 
to  the  sound  source. 
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5-   Signal  Data  Recorder-Reproducer  Model  AN/lMQ-8(v).   The 
AN/UNQ-8(v)  utilized  in  this  experimental  work  was  modified  for  tape 
loop  operation.   An  external  tape  loop  transport  was  added  and  the 
record  circuitry  was  modified  in  order  to  disable  the  normal  record  push 
button  and  to  provide  for  a  remote,  hand-held  record  switch.   In  the  work 
reported  herein,  the  recorder-reproducer  was  operated  in  direct  record  in 
the  loop  mode  at  a  tape  speed  of  30  in/sec.   The  recorder-reproducer 
channel  utilized  was  calibrated  as  directed  in  the  AN/UNQ-8(v)  manual 
prior  to  its  use.   The  3  dE-  down  points  of  the  frequency  response  of  the 
AN/UNQ-8(V)  used  as  indicated  above  are  50  Hz  and  60  kHz. 

6.   ELGENCO  Random  Noise  Generator  Model  603A.   We  investigated 
the  instantaneous  amplitude  distribution  of  the  output  of  the  Model  603A. 
The  Model  603A  was  to  be  the  source  of  our  noise  waveforms  and  we  wanted 
to  establish  whether  the  source  waveform  were  truly  random. 

The  output  of  the  Model  603A  was  analyzed  using  a  standard 
statistical  program  called  SEQSTAT  available  on  the  computer  library  at 
the  Applied  Research  Laboratories.   The  output  was  alias  filtered  at 
10  kHz  using  a  low  pass  filter  of  slope  48  dB/octave.   The  filter  output 
was  then  sampled  at  25.22  kHz  using  an  analog-to-digital  converter.   Four 
sample  sequences  50,000  data  points  in  length  were  formed.   In  SEQSTAT 
sequences  of  5000  data  samples  are  considered  and  the  mean,  variance, 
skew,  and  kurtosis  of  the  sample  record  are  computed.   In  addition  a 
histogram  is  formed  and  compared  to  a  gauss ianly  distributed  variable 
having  mean  and  variance  equal  to  those  calculated  from  the  data. 


All  external  and  internal  modifications  to  the  AN/UNQ-8(v)  had  been 
designed  and  built  by  R.  K.  Goodnow  of  the  Applied  Research  Laboratories 
staff. 
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Consider  a  random  variable  x  having  a  probability  distribution  function 
f(x).   The  definitions  of  the  first  four  central  moments  are  [see,  for 
example,  Wilks  (1963)] 

Mean         m  =  E(x),  (5-1 ) 

2  2 

Variance     o     =   E[ (x-m)  ],  (5-2) 

Skew        y±   =  E[ (x-m)5]/a5,  and  (5-5) 

Kurtosis     70  =  E[  (x-m)  ]/a  ,  (5-^) 

where  the  notation  E  means  expected  value.   For  a  gaussianly  distributed 
random  variable  the  theoretical  value  of  7  and  7  are 

71  =  0 

r2  =  3-0   .  (5-5) 

The  first  four  moments  of  the  Model  603A  output  were  calculated  for  12 
samples  records  of  5000  samples  each  and  the  moments  of  each  sample 
record  were  then  ensemble  averaged.   The  results  were  as  follows: 

m  =  0.0039  V, 

a  =  1.227  V, 

(5-6) 
7  =  0.012,  and 

72  =  2.92. 

The  true  rms  voltage  at  the  output  of  the  alias  filter  was  1.22  V, 
agreeing  very  well  with  o.      The  ensemble  averaged  histogram  of  the 
sampled  data  agreed  closely  with  the  gaussian  distribution  curve  based 
on  (m, c)  in  result  Eq.  (5-6). 

We,  therefore,  conclude  that  the  ELGENCO  Model  603A  produced 
nearly  gaussianly  distributed  noise  when  operated  in  the  5  Hz  to  10,000  Hz 
range . 
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7.  DuKane  Power  Amplifier  Medallion  Model  1A921.   This 
device  is  power  rated  by  the  manufacturer  at  200  W  true  rms  and  500  W 
peak  instantaneous  power.   The  Model  1A921  -was  used  to  power  the 
ID-75  Horn  Driver.   The  maximum  usable  output  from  the  Model  1A921 
■was  set  by  the  capacity  of  the  ID-75  ("when  operated  in  the 

pulsed  mode).   The  additional  power  capability  of  the  Model  1A921  was 
used  to  minimize  clipping  of  the  amplifier  output.   If  we  assume  that  the 
amplifier  input  is  gaussianly  distributed  random  noise  and  the  true  rms 
amplifier  output  is  75  W,  then  the  probability  of  waveform  clipping, 
obtained  from  the  gauss ian  distribution  function  tables  [see,  for  example, 
Baird  (196A)],  is  approximately  \°jo.      The  3  dB  down  points  of  the  Model 
1A921  frequency  response  are  at  20  Hz  and  20  kHz.   Hence,  the  Model  1A921 
did  not  limit  the  source  bandwidth  used  in  the  transmit  systems  in  our 
noise  experiments.   In  fact,  as  has  already  been  mentioned,  the  Krohn-Hite 
filter  was  used  to  limit  the  signal  to  the  amplifier  to  the  bandwidth 
usable  by  the  horn  driver.   The  rated  total  harmonic  distortion  of  the 
Model  1A921  is  less  than  1  l/2$>  when  operated  at  output  powers  less  than 
200  W  rms.   Since  we  operated  at  powers  far  less  than  200  W,  it  is 
unlikely  that  the  Model  1A921  output  was  significantly  distorted.   In 
addition,  we  monitored  the  amplifier  output  using  an  oscilloscope  during 
our  experiments  and  never  observed  waveform  distortion.   Since  we  are 
interested  in  the  change  in  the  waveform, the  presence  of  distortion, 
had  any  been  present,  would  not  have  concerned  us. 

8.  University  Horn  Driver  ID-75.   The  ID-75  has  a  bandpass 
filter  frequency  characteristic.   When  the  ID-75  is  operated  at  75  W 
output  the  3  dB  down  points  of  the  passband  are  approximately  700  Hz 
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and  2200  Hz.   The  ID-75  frequency  characteristic  essentially  sets  the 
bandwidth  at  the  source  for  our  noise  experiments. 

The  ID-75  is  rated  for  75  W  rms  continuous  output.   All  of  our 
high  intensity  (approximately  160  dB)  results  were  obtained  while  pulsing 
the  transmit  system.   We  found  that  higher  overall  source  SPL's  could 
be  reached  in  the  pulse  mode.   We  estimate  that  the  ID-75  in  the  pulse 
mode  can  safely  (i.e.,  nondestructively)  accommodate  75  W  rms  electrical 
power.  At  this  power  level,  however,  the  acoustic  output  waveforms  are 
somewhat  distorted.   Examples  of  the  distortion  caused  by  the  driver  will 
be  presented  in  Chapter  6.   The  amount  of  horn  driver  induced  distortion 
does  not  concern  us  here  because,  as  we  mentioned  in  the  discussion  of 
the  power  amplifier,  we  are  interested  in  the  change  in  the  waveform  with 
propagation.   Hence,  all  we  require  at  the  source  is  that  the  waveform, 
whatever  it  may  be,  be  known. 

We  investigated  the  output  of  the  ELGENCO  noise  generator  and 
established  that  its  output  was  gaussian.   Downstream  of  the  noise  genera- 
tor, however,  we  have  significantly  altered  the  signal  by  filtering  and 
amplifying  it.   The  filtering  and  amplifying  action  are  linear  trans- 
formations and,  as  is  well-known,  gaussian  processes  are  invariant  under 
a  linear  transformation  [see,  for  example,  Papoulis  (1965)].   Hence  in 
the  absence  of  nonlinear  distortion  the  output  of  the  horn  driver  is  also 
gaussian.   We  know,  however,  that  the  horn  driver  at  high  output  levels 
does  nonlinearly  distort  the  signal  somewhat.   Hence  the  distribution 
function  of  the  horn  driver  output  may  not  be  exactly  gaussian.   The 
only  effect  that  this  possible  change  in  the  distribution  function  has 
in  regard  to  our  computations  is  its  effect  on  our  calculations  of  spectra. 
The  methods  used  depend  on  the  noise  amplitudes  being  normally  distributed. 
We  acknowledge,  therefore,  that  since  the  horn  driver  at  high  power  does 
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nonlinearly  distort  our  signal  somewhat  the  calculation  of  spectra  may 
be  affected.   However,  we  feel  confident  that  the  effect  is  of  a  very 
minor  nature. 

9.   General  Radio  Pulse  Generator  Type  1396-A.   The  Type  I396-A 
was  used  in  all  the  finite-amplitude  propagation  experiments  to  pulse  the 
transmit  system.   It  was  not  used  in  the  measurement  of  the  pipe  termina- 
tion reflection  coefficient  "because  during  the  "off"  portion  of  its  duty 
cycle  the  Type  I396-A  has  a  feedthrough  only  40  dB  below  the  output  during 
the  "on"  portion  of  the  cycle.   Since  we  were  interested  in  reflections 
more  than  40  dB  below  the  incident  signal,  the  Type  I396-A  was  unsatis- 
factory for  use  in  reflection  measurements. 

10.  Wavetek  Function  Generator  Model  144.   This  device  was  used 
to  pulse  the  transmit  system  during  the  reflection  coefficient  measurements. 
The  Model  144  exhibits  zero  feedthrough  during  the  "off"  portion  of  its 
duty  cycle. 

11.  Daven  Attenuator  Type  T64l.   This  device  was  used  at  the 
Model  1A921  power  amplifier  input  to  control  the  amplifier  output  level. 

12.  United  Variable  Tap  Transformer.   The  transformer  was  used 
in  a  1:1  configuration  to  isolate  the  input  of  the  power  amplifier  from 
the  laboratory  electrical  ground. 

b.   The  Pipe 

The  pipe  consists  of  eight  12  ft  sections,  seven  flanges,  one 
driver  mount  and  a  nonreflecting  termination.   The  pipe,  less  the  termina- 
tion, was  originally  constructed  for  the  experiments  reported  by  McKittick, 
Blackstock,  and  Wright  (1967);  certain  modifications  have  been  made  for  the 
present  work.   The  termination  is  constructed  of  fiberglass  and  the  pipe 
sections  and  fittings  of  6061-T6  AL  aluminum. 
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1.  The  Pipe  Sections.   The  pipe  section  internal  diameter  is 

2  in.  and  wall  thickness  is  1/8  in.   The  ends  of  each  section  were  squared 
and  the  outside  surfaces  within  3  in.  of  the  ends  were  turned  round  on  a 
lathe  to  facilitate  a  close  fit  with  the  flanges. 

2.  The  Pipe  Flanges.   An  assembly  drawing  of  one  of  the  pipe 
flanges  appears  in  Fig.  (5-2).   The  flanges  were  designed  to  minimize 
reflections  at  the  pipe  joints.   Each  flange  has  a  probe  hole  sized  to 
accept  a  l/k   in.  microphone.   The  microphone  body  is  mounted  in  a  teflon 
holder  that  seats  on  the  shoulder  in  the  microphone  hole  [see  Fig.  (5-2)]. 
Since  the  flanges  are  not  of  equal  diameter,  the  depth  of  this  shoulder 
varies  in  such  a  fashion  that  the  face  of  the  mounted  microphone  is  always 
flush  with  the  inside  surface  of  the  pipe.   The  teflon  holder  also 
reduces  the  mechanical  pickup  of  signals  propagating  in  the  pipe  walls. 
When  not  in  use,  the  microphone  hole  is  blocked  by  an  aluminum  plug  whose 
end  is  flush  with  the  inside  surface  of  the  pipe.   The  taper  pins  ensure 
exact  alignment  of  the  flanges  each  time  the  flanges  are  bolted  together. 
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FIGURE  5  2 
AN  ASSEMBLY  DRAWNG  OF   A  PIPE  FLANGE 
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3.   Horn  Driver  Mount.   The  horn  driver  mount  consists  of  a 
flange  similar  to  that  shown  in  Fig.  (5-2).   On  the  driver  side  of  the 
flange  an  8  in.  conical  transition  piece  is  used  to  change  the  internal 
diameter  from  2  in.  to  1  in.  to  match  the  throat  diameter  of  the 
University  series  horn  drivers. 

k.        The  Temination.   The  pipe  termination  was  built  as 
suggested  in  the  paper  by  Burns  (1971).   Of  Burns'  three  termination 
designs  we  elected  to  use  the  full  wave  tangent  taper.  The  terminat- 
tion  consists  of  two  sections  which  are  installed  in  series  in  the  final 
length  of  pipe  section.  The  first  section  of  the  termination  encountered 
by  the  propagating  wave  consists  of  a  gradually  tapered,  one  meter  long 
wedge.  The  second  section,  also  one  meter  long,  completely  fills  the 
pipe  cross  section.  The  density  of  the  fiberglass  used  is  O.O503  gm/cm  . 
The  leading  edge  of  the  termination  is  located  approximately  5  ft  down- 
stream of  the  seventh  (and  last)  flange.  Therefore,  the  termination 
did  not  interfere  with  readings  taken  at  the  seventh  flange. 

The  reflection  coefficient  of  the  pipe  termination  was 
measured  by  the  pulse  method.   In  Table  I  we  present  the  reflection 
coefficient,  corrected  for  the  round  trip  loss  due  to  tube  wall 
attenuation  between  the  measurement  point  and  the  termination. 
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TABLE  I 
Reflection  Coefficients  for  the  Pipe  Termination 

Frequency  (Hz)  Reflection  Coefficient  (dB) 

500  -39 . 8 

750  -35.7 

iooo  -37.2 

1500  -I+0.8 

2000  -42.1 

2500  -In.  7 

3000  -^1.8 

3500  -1+2.0 

1+000  -i+i.  9 

Burns  (1971)  notes  that  a  reflection  coefficient  appreciably  below 
-1+0  dB  is  hard  to  achieve  in  practice, 
c.   The  Receiving  System 

The  accuracy  of  the  receiving  system  calibration  established 
the  accuracy  of  our  results.   Each  system  component  was  calibrated. 
The  accuracy  of  measurements  made  with  the  various  devices  was  determined 
during  the  device  calibration.   In  this  subsection  "we  report  the  measure- 
ment standard  deviations  of  the  various  components  of  the  system. 
Microphone  sensitivity  H  "was  measured.   The  values  of  H  reported  here 
are  for  the  microphone  assemblies,  i.e.,  the  microphone,  preamplifier, 
and  cable,  connected  as  they  were  used  in  our  experiments. 

1.  Bruel  and  Kjaer  Preamplifiers  Type  2619  and  Type  2615. 
The  Type  26l5  "was  used  only  in  the  microphone  calibration.   The  Type  2619 
preamplifiers  were  used  for  all  data  collection. 

2.  Bruel  and  Kjaer  Cable  Type  A00029.   This  cable  is  30  m 
in  length. 
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3.   Laboratory  Constructed  Microphone  Cable.   This  device 
was  used  in  conjunction  with  the  source  microphone.  The  cable  type  is 
RG-62  B/U.   The  end  fittings  are  of  Bruel  and  Kjaer  manufacture.  A 
careful  check  was  made  to  ensure  that  signals  detected  by  the  source 
microphone  were  not  degraded  by  use  of  this  cable.  At  frequencies  up  to 
70  kHz  no  signal  distortion  or  loss  of  amplitude  was  observed.   The 
cable  is  100  ft  in  length. 

k.        Bruel  and  Kjaer  l/k   in.  Microphones  Type  ^136.   The  \/h   in. 
microphones  were  calibrated  by  means  of  a  comparison  method.   Three 
Bruel  and  Kjaer  1/2  in.  microphones  were  calibrated  by  the  chamber 
reciprocity  method  [see  Bruel  and  Kjaer  (1965)].   The  1/2  in.  microphones 
were  used  to  calibrate  a  Bruel  and  Kjaer  Type  4220  pistonphone.   The 
pistonphone  was  then  used  to  calibrate  the  l/h   in.  microphones.   In 
Table  II  we  present  the  measured  l/k-   in.  microphone  assembly  sensitivi- 
ties. A  microphone  assembly  consists  of  a  microphone,  a  preamplifier,  and 
a  cable  connected  in  series.   The  values  reported  are  the  sensitivities 
at  250  Hz,  76O  mm  Hg,  and  200  Vdc  polarization  voltage. 

TABLE  II 
Microphone  Assembly  Sensitivities 

H  (dB  re  1  V/|ibar) 

source  assembly  -77*2 

remote  assembly  -76.8 

In  the  experiments  reported  in  this  work  \/h   in.  microphone  serial  I76OO5, 
Type  2619  preamplifier  serial  32725I+,  and  the  locally  constructed  cable 
were  used  to  monitor  the  sound  at  the  source.   Microphone  serial  I76OO7, 
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Type  2619  preamplifier  serial  386559,  and  Cable  A00029  were  used  at 
the  remote  locations.  After  the  completion  of  the  initial  calibration, 
the  microphone  assembly  sensitivities,  checked  using  the  calibrated 
pistonphone,  were  repeatable  within  ±0.15  dB.   The  nominal  frequency 
response  characteristic  of  the  l/k   in.  microphones  is  flat  within  ±1  dB 
over  the  range  of  50  Hz  to  70  kHz.   Series  addition  of  the  preampli- 
fier and  cable  for  both  the  remote  and  source  microphone  assembly 
had  no  observable  effect  on  the  microphone  frequency  response.   The 
time  domain  characteristics  of  the  l/k   in.  microphones  made  them  well 
suited  to  meet  our  experimental  needs.   The  nominal  rise  times  are  very 
short  (approximately  6  usee  in  response  to  a  step  input)  and  the  micro- 
phones exhibit  practically  no  ringing  when  installed  as  they  were  in  our 
experiments . 

Use  of  a  finite  size  microphone  requires  that  we  investigate 
the  integration  effect  associated  with  measuring  a  sound  field  at  a 
grazing  incidence.   Mounted  as  used  in  the  pipe,  the  microphone  is  simply 
a  baffled  circular  piston  receiver.   The  directivity  characteristic  D(k, 0) 
is  easily  calculated  to  be  [see,  for  example,  Kinsler  and  Frey  (1950)] 

2J  (ka  sin  Q) 

where  a  =  the  microphone  radius  and  9   is  the  angle  of  incidence  with 
respect  to  the  normal  to  the  microphone  face.   In  the  present  case 
0=90°  and  Eq.  (5-l)  reduce  to 


2 J  (ka) 
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For  the  following  reasons  we  did  not  correct  our  reported  theoretical 
waveforms  for  microphone  directivity: 

1)  The  3  cLB  down  point  of  D(k)  for  a  \jh  in.  microphone  is 
27.5  kHz.  Only  a  very  small  portion  of  the  energy  in  our  measured  and 
computed  spectra  occurs  at  frequencies  above  27.5  kHz. 

2)  Test  runs,  designed  to  discover  the  effect  of  the 
inclusion  of  D(k),  were  made.   In  these  tests  the  computation  algorithm 
was  modified  to  include  D(k)  in  the  calculation  of  theoretical  waveforms. 
The  corrected  waveforms  (i.e.,  corrected  for  microphone  directivity)  were 
then  compared  to  uncorrected  waveforms  computed  from  similar  source 
data.   The  comparison  showed  that  the  only  difference  between  the  two  sets 
of  data  was  a  small  (5  to  6  M-sec)  increase  in  shock  rise  time  in  the 
corrected  data.   The  difference  in  shock  rise  time  is  less  than  the 
experimental  uncertainty  in  the  recording  of  our  data. 

5-   Bruel  and  Kjaer  Selector  Switch  Type  4^08.  Use  of  this 
device  allowed  us  to  easily  switch  back  and  forth  between  the  two 
microphones.  Use  of  the  2-channel  selector  switch  had  no  adverse 
effect  upon  signal  measurement. 

6.  Bruel  and  Kjaer  Power  Supply  Type  2803. 

7.  Bruel  and  Kjaer  Pistonphone  Type  4220.   This  device  was 
used  to  check  calibrate  the  microphone  assemblies  before  every  data 
collection. 

8.  Tektronix  Oscilloscope  Model  5^5B  with  Plug-in  Type  1A1. 
The  measurement  standard  deviation  of  the  model  5^5B  oscilloscope  after 
calibration  was  ±0.25  dB.   The  Model  5^5B  was  used  for  several  signal 
monitoring  purposes  that  are  outlined  briefly  in  Section  2. 
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9.   General  Radio  Wave  Analyzer  Type  19OO-A.   This  device 
■was  used  in  the  collection  of  noise  wave  spectra.   For  all  spectral 
measurements  reported  in  this  work  the  50  Hz  crystal  was  used.   The 
nominal  measurement  standard  deviation  of  the  Type  1900-A  used  with  the 
Recorder  Type  1521B  was  ±0.2  dB. 

10.  General  Radio  Recorder  Type  1521B.   This  pen  type 
recorder  was  used  in  conjunction  with  the  wave  analyzer  to  record  spectra. 

11.  Hewlett  Packard  Oscilloscope  Model  18IAIAR  with  ^-Channel 
Vertical  Amplifier  Type  iQoh-A   and  Time  Base  and  Delay  Generator 

Model  l821A.   This  memory-type  oscilloscope  was  used  in  the  noise  experi- 
ments.  Oscillograms  of  the  noise  pulses  were  made  using  the  memory 
feature  of  the  Model  18IAIAR.  At  the  time  that  this  device  was  used, 
it  was  new  and  had  recently  been  factory  calibrated.   The  observed 
measurement  standard  deviation  of  the  Model  18IAIAR  was  ±0.1  dB. 

12.  Bruel  and  Kjaer  Filter  Type  l6l0.   This  l/3  octave  band 
filter  set  was  used  in  the  measurement  of  small-signal  attenuation 
coefficients . 

13.  Bruel  and  Kjaer  Microphone  Amplifier  Type  260^.   This 
amplifier  was  used  only  in  the  measurement  of  small-signal  attenuation 
coefficients .  We  measured  the  standard  deviation  of  the  response  of 
the  Type  26oh   and  obtained  ±0.3  dB. 

3.   Experimental  Procedures 

We  conducted  three  types  of  experiments  in  the  research 
reported  here .   The  experiments  were  as  follows : 

a.  Measurement  of  small-signal  attenuation  coefficients. 

b.  Measurement  of  deterministic  waveforms. 

c.  Measurement  of  noise  waveforms. 
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We  shall  "briefly  discuss  our  procedure  and  identify  the  equipment  used 
in  each  of  the  experiments.   The  equipment  -will  be  identified  by  its 
system  number,  e.g.,  the  General  Eadio  Oscillator  Type  lJlO-A  -will  be 
referred  to  as  item  a-1.   Since  the  pipe  -was  used  in  all  three  experi- 
mental procedures  it  will  not  be  listed  in  the  individual  experimental 
descriptions.   Before  each  experiment  the  microphone  assembly  sensitivities 
were  checked  using  the  pistonphone. 

a.  Measurement  of  Small-Signal  Attenuation  Coefficients. 
The  devices  used  in  this  experiment  were  items  a-2,  a-3,  a-^,  a-7,  a-8, 
a-11,  a-12,  c-1,  c-2,  c-3,  c-*+,  c-5,  c-7,  c-8,  c-12,  c-lj.   The  purpose 
of  this  experiment  was  to  provide  attenuation  coefficients  for  use  in 
the  computation  algorithm.   The  procedure  used  consisted  of  generating 
a  low  amplitude  (typically  105  dB)  sinusoid  at  the  source  and  recording 
the  signal  level  at  the  various  remote  locations.   The  B  and  K  oscillator 
and  the  Beckman  timer  were  used  to  obtain  maximum  frequency  stability 
and  accurate  frequency  measurement,  respectively.   The  oscilloscope 

was  used  to  monitor  source  electrical  and  acoustic  signals  and  acoustic 
signals  at  the  remote  locations.  The  attenuation  coefficients  were 
calculated  by  plotting  the  observed  sound  pressure  levels  as  a  function 
of  distance  and  obtaining  the  best  fit  line  by  linear  regression. 

b.  Measurement  of  Deterministic  Waveforms. 

The  items  used  were  a-1,  a-3,  a-U,  a-7,  a-8,  a-9,  a-11,  a-12,  c-1,  c-2, 
c-3,  c-4,  c-5,  c-6,  c-7,  and  c-8.  The  basic  procedure  used  consisted 
of  pulsing  the  source  at  a  high  amplitude  and  observing  the  waveform 
at  the  various  pipe  locations.  At  each  location  an  oscillogram  was 
made  using  the  Model  5^5B  oscilloscope  camera.   The  General  Radio 
Oscillator  was  used  instead  of  the  B  and  K  oscillator  because  the  former 
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is  compatible  with  the  General  Radio  Pulse  Generator.   The  Model  5^5B 
oscilloscope  was  also  used  to  monitor  the  electrical  signal  to  the  horn 
driver  and  the  source  acoustic  signal  in  order  to  ensure  that  the  source 
signal  remained  stable  during  a  particular  data  collection  period. 

c.   Measurement  of  Noise  Waveforms.   The  items  used  -were 
a-1,  a-^,  a-5,  a-6,  a-7,  a-8,  a-9,  a-11,  a-12,  c-1,  c-2,  c-3,  c-k, 
c-5,  c-6,  c-7,    c-9,  c-10,  c-11.   This  experiment  is  somewhat  more 
complicated  than  experiments  (a)  and  (b). 

In  the  noise  -waveform  experiments  a  pulse  of  filtered 
electrical  noise  was  recorded  on  an  endless  tape  loop.   The  filter 
bandpass  was  500  "to  3500  kHz.   The  recording  was  made  by  gating  the 
output  of  the  ELGENCO  noise  generator  by  use  of  the  General  Radio  Pulse 
Generator  and  General  Radio  Oscillator.   The  pulse  was  recorded  using 
the  AN/UNQ-8( V)  modified,  as  mentioned  earlier,  for  tape  loop  operations. 
The  AN/UNQ,-8( V)  was  then  used  to  reproduce  the  noise  pulse.   The  pulse 
was  amplified  and  applied  to  the  ID- 75  Horn  Driver.   The  acoustic  noise 
pulse  was  then  picked  up  by  the  remote  microphone  and  the  signal 
recorded  by  the  Model  181A1AR  oscilloscope.   This  oscilloscope  was 
also  used  to  monitor  and  record  the  source  signal. 

In  the  noise  spectra  experiments  the  ELGENCO  Model  603A 
output  was  bandpass  filtered,  amplified,  and  applied  continuously  (not 
pulsed)  to  the  ID-75.   Again,  the  filter  bandpass  was  500  to  3500  kHz. 
The  maximum  sound  pressure  level  that  could  be  safely  obtained  in  the 
continuous  wave  (cw)  mode  was  approximately  15^  dB.   The  General  Radio 
Wave  Analyzer  and  General  Radio  Recorder  were  used  to  record  spectra 
at  the  various  pipe  measurement  stations. 
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k.        Experimental  Limitations 

Chapter  3  concluded  -with  a  listing  of  several  assumptions  on 
which  the  formulation  of  modified  •weak-shock  theory  is  based.  For  the 
convenience  of  the  reader,  we  shall  restate  these  assumptions: 

a.  Only  plane  waves  are  present. 

b.  The  wave  amplitude  is  appropriate  for  the  use  of  modified 
weak-shock  theory. 

c.  The  pipe  is  a  wide  pipe,  i.e., 


A  «  R  «  -|- 
k  A 


d.  The  pipe  is  circular,  infinitely  long,  and  absolutely  rigid, 

e.  The  enclosed  fluid  medium  is  a  homogeneous  perfect  gas. 
In  this  section  we  shall  explain  the  physical  limitations  placed  on  our 
experiments  by  the  assumptions  stated  above. 

The  plane  wave  assumption  establishes  an  important  operating 
restriction.   The  cutoff  frequencies  of  nonplanar  modes  in  a  circular 
pipe  are  inversely  proportional  to  the  pipe  radius.   Hence,  the  selection 
of  the  pipe  i.d.  establishes  for  practical  purposes  the  upper  limit  on 
the  source  frequency.   The  cutoff  frequencies  are  given  by  the  relation 

P  c 


'»-lr     •  (5"5) 


In  Eq.  (5-3)  P   is  the  nth  root  of  the  equation 


mn 


j  ■  (PJ  =  0  (J-*) 
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■where  J  '  =  the  derivative  of  the  nth  order  Bessel  function.   The  first 
several  cutoff  frequencies  are  given  in  Table  III.   The  row  and  column 
indices  in  Table  III  identify  the  root  of  Eq.  (5-4)  used  in  Eq.  (5-5) 
to  calculate  the  given  frequency. 

TABLE  III 
Cutoff  Frequencies  (Hz)  in  a  2  in.  i.d. 
Circular  Pipe 


m  =  0 

m  =  1 

m  =  2 

n  =  0 

0 

8280 

15170 

n  =  1 

394o 

11510 

18410 

n  =  2 

6590 

1^500 

21550 

The  (0,0)  mode  is  planar.   The  (0,  l)  and  (0,2)  modes  are  spinning 
modes  and  the  (l,0)  mode  is  radial.   In  order  to  avoid  the  generation 
of  nonplanar  modes  at  the  source,  source  frequencies  in  the  experiments 
reported  here  were  maintained  below  39^0  Hz. 

The  waveform  steepening  caused  by  nonlinear  propagation  shifts 
energy  to  frequencies  well  above  the  39^0  Hz  limit.   Does  this  energy 
transferred  to  high  frequencies  generate  nonplanar  modes?  We  are  unable 
to  conclusively  answer  this  question,  but  we  feel  that  it  does  not.   Our 
experimental  data  and  that  of  other  researchers  [see,  for  example,  Pernet 
and  Payne  (1969)]  indicates  that  high  frequency  components  generated  by 
nonlinear  steepening  of  plane  waves  remain  planar.   The  planar  shock  seems 
to  represent  a  very  stable  mode  of  propagation  in  a  tube. 

We  have  assumed  a  plane  particle  velocity  wave  whereas  in 
reality  the  particle  velocity  must  vanish  at  the  pipe  wall.   Hence,  the 
particle  velocity  phase  front  is  plane  only  outside  the  boundary  layer. 
We  measure  pressure,  however,  and  the  pressure  is  planar  across  the  pipe 
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cross  section  through  the  boundary  layer.   Therefore,  measurements  made 
with  microphones  flush  with  the  tube  wall  are  completely  acceptable. 

Modified  weak-shock  theory  is  applicable  if  the  waves  being 
considered  are  not  too  strong.   Blackstock  (1964)  maintains  that  waves 
having  sound  pressure  levels  less  than  174  dB  may  be  treated  by  weak- 
shock  theory.   Temkin  (1969),  in  a  more  detailed  analysis,  limits  weak- 
shock  theory  to  165  dB.   In  our  work  the  maximum  sound  pressure  level 
reached  was  162  dB.   In  Appendix  E  we  present  our  analysis  of  the  maximum 
amplitude  that  may  be  considered  using  weak-shock  theory.   Our  results 
indicate  that  165  dB  is  a  realistic  and  practical  figure. 

The  wide  pipe  assumption  does  not  impose  any  restrictions  on 
our  experiments  over  and  above  those  restrictions  already  imposed  by  the 
limitations  of  our  equipment.   The  low  frequency  limit,  set  by  the  ID-75 
response,  is  approximately  700  Hz.   The  high  frequency  limit,  set 
by  the  microphone  frequency  response,  is  27.5  kHz.   Within  the  frequency 
range  of  700  Hz  to  well  above  27.5  kHz  our  pipe  is  a  wide  pipe. 

The  upper  limit  of  the  wide  pipe  assumption,  i.e.,  R«Tl/(k  A), 
was  derived  by  comparing  Kirchhoff  theory  to  thermoviscous  losses  in  the 
pipe  mainstream.   We  ignored  other  loss  mechanisms.   If  we  include  relaxa- 
tion in  the  mainstream  attenuation  and  compare  the  boundary  layer  attenuation 
coefficient  and  the  mainstream  attenuation  coefficient,  we  find  that  the 
frequency  at  which  the  two  coefficients  approach  each  other  is  significantly 
lower  than  in  the  case  where  relaxation  effects  are  ignored.   Using  the 
results  of  Bass,  Evans,  and  Sutherland  (1970),  we  have  made  the  above 
comparison  and  found  that  for  frequencies  less  than  100  kHz  the 
boundary  layer  attenuation  coefficient  exceeds  the  mainstream  term. 
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Hence,  we  conclude  that  we  incur  negligible  error  in  our  results  by 
ignoring  relaxation  effects. 

The  final  assumption  is  a  consequence  of  our  use  of  Kirchhoff 
theory  and  also  of  the  simple  wave  flow  assumption.   The  pipe  termination 
is  meant  to  make  the  pipe  free  of  reflection  and,  in  this  sense,  simulate 
an  infinitely  long  pipe.   The  figures  in  Table  I  make  clear  that  the  pipe 
is  indeed  a  progressive  wave  tube. 
5.   Analysis  of  the  Expected  Measurement  Accuracy 

The  accuracy  of  the  absolute  sound  pressure  level  measurements 
reported  in  this  work  may  be  easily  estimated.   If  we  assume  that  the 
various  errors  are  independent  and  that  each  is  normally  distributed, 
the  standard  deviation  of  the  measurement  is 


S  =  [a2  +  b2  +  c2  +  ....]1/2    ,  (5-5) 


where  a,  b,  c,  etc.,  are  the  standard  deviations  of  responses  of  the  varioi 
components  used  to  make  a  particular  measurement  and  S  is  the  standard 
deviation  of  the  overall  system  response.   The  measurement  standard  devia- 
tions for  the  equipment  used  in  each  experiment  are  listed  below.   Values 
marked  by  the  symbol*  are  worst  case  values  obtained  from  the  various 
equipment  technical  manuals.   All  other  values,  not  so  marked,  were 
measured.   The  measurement  standard  deviations  are 

a)  Microphone  assembly  calibration:   O.15  dB. 

b)  Pistonphone  output:   0.2  dB* 

c)  Voltage  calibration  of  the  microphone  amplifier:   0.2  dB. 

d)  Reading  uncertainty:   0.1  dB. 

e)  Polarization  voltage  deviation  from  200  V:   0.2  dB*. 

f )  Uncertainty  in  the  frequency  response  of  the  measuring 
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equipment:   0.3  dB  . 
g)  Voltage  calibration  of  the  Tektronix  Oscilloscope  5^5B: 

0.25  dB. 
h)  Voltage  calibration  of  the  HP  181A1AR  oscilloscope:   0.1  dB. 
i)  Voltage  calibration  of  the  General  Radio  Wave  Analyzer  and 
Recorder:   0.2  dB  . 
In  Section  2  we  listed  the  equipment  used  in  each  experiment.   Using 
Eq..  (5-5)  and  the  standard  deviations  listed  above,  we  may  calculate  the 
measurement  standard  deviation  for  the  three  experiments  as  follows: 

a.  Measurement  of  Attenuation  Coefficients.   The  measurement 

standard  deviation  S  is 

a 

r  2   ,  2    2    n2    2   J2.-X/2 
S  =  [a  +  b  +c  +d  +  e  +  f  ]  ' 

(5-6  ) 
=  ±0.54  dB 

b.  Measurement  of  Deterministic  Waveforms.   The  measurement 

standard  deviation  S.  is 

b 


r  2   ^2   A2          2    _2  2-.1/2 

S  =  [a  +B  +d  +e  +f  +g]7 

=  ±0.55  dB 

c.   Measurement  of  Noise  Waveforms.  The  measurement 


(5-7) 


standard  deviation  S  is 

c 


S 
c 


r  2   J2.        ,2    2   ^2   ,.2.1/2 
=  [a  +b  +d  +e  +  f  +h  J  ' 


(5-8) 


=  ±0.5  dB 
Noise  spectra  were  also  measured.   The  standard  deviation  for 
the  spectral  measurements  is  obtained  by  replacing  item  h  in  Eq.  (5-8) 
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by  item  j  from  the  list  of  standard  deviations.   The  resulting  value  of 
S  is  ±0.53  cLB. 

Since  the  three  accuracy  numbers  are  all  about  the  same,  we  can 
safely  say  that  the  standard  deviation  of  the  measurements  is  ±0.5  dB. 
If  we  refer  to  the  cumulative  distribution  function  for  a  gaussian  random 
variable  [see,  for  example,  Baird  (1964)],  we  see  that  the  value  of  S  may 
be  interpreted  as  measuring  that  there  is  a  68fo  chance  of  obtaining 
errors  less  than  S,  a  90/0  chance  of  obtaining  errors  less  than  1. 65  S,  and 
a  99fo  chance  of  obtaining  errors  less  than  2.58  S.   Thus,  with  a  standard 
error  of  ±0.5  dB  the  probability  of  obtaining  errors  larger  than  ±0.8  dB 
is  about 


CHAPTER  6    EXPERIMENTAL  AND  THEORETICAL  RESULTS 

1.   Introduction 

In  Chapter  6  we  present  and  discuss  the  theoretical  and 
experimental  results  that  we  have  obtained.   The  chapter  is  divided 
into  several  subsections,  each  devoted  to  a  particular  group  of  results. 
The  subsections  are  as  follows. 

1.  Introduction. 

2.  Measurement  of  tube  wall  attenuation  coefficients. 

3.  Results  concerned  with  waves  that  are  sinusoidal  at 
the  source. 

h.      A  comparison  of  results  obtained  using  waves  that  are 
sinusoidal  at  the  source  with  those  of  other  researchers. 

5.  Noise  waveform  results . 

6.  Predictions  and  measurement  of  noise  spectra. 

7-   The  relation  of  our  noise  waveform  results  to  results 
obtained  using  the  "Burgers  turbulence"  model. 

Prior  to  each  experimental  measurement,  the  environmental 

conditions  in  the  laboratory  were  carefully  measured.   The  laboratory 

building  is  maintained  at  72°  ±2°F.   The  realtive  humidity,  although 

not  specifically  controlled,  also  tends  to  be  quite  constant.   The  range 

of  relative  humidities  observed  during  the  course  of  our  experiments 

was  from  40$>  to  60</0.   Building  pressure  varies  with  outside  pressure. 

Nearly  constant  temperature  means  that  the  small-signal  sound  speed 

c  is  nearly  constant. 
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The  experimental  results  described  here  were  obtained  by 
making  measurements  at  the  eight  pipe  probe  holes.   We  refer  to  the 
distances  from  the  source  to  these  locations  using  the  nearest  integer 
number  of  feet.   The  exact  distance  from  the  horn  driver  diaphragm  to 
the  various  measurement  locations  is  given  in  Table  IV. 

TABLE  IV 

Measurement  Locations 

Measurement  Exact  Distance 

Location  from  Source 
(ft) 

1  0.855 

2  13.015 

3  25.065 
A  37-175 

5  ^9.255 

6  61.355 

7  73.^35 

8  85.495 

In  the  theoretical  results  the  distances  indicated  are  the 
exact  distances  used  in  the  computations.   The  very  snail  differences 
between  the  propagation  distances  used  in  the  theoretical  calculations 
and  the  actual  distances  shown  in  Table  IV  are  insignificant.   The  wave- 
form at  the  point  x=l  ft  is  referred  to  as  the  source  waveform  in  the 
results  discussed  below.   Although  the  true  source  is  the  horn  driver 
diaphragm,  we  are  not  concerned  here  with  the  actual  waveform  at  the 
diaphragm.   Rather  we  are  concerned  with  the  change  in  the  waveform  as 
it  propagates.   Hence,  for  our  purposes  we  may  regard  the  1  ft  measure- 
ment location  as  the  source  point.   The  use  of  a  stationary  source 
location  is  consistent  with  assumptions  made  in  the  derivation  of 
Eq.  (3-6). 
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In  several  of  the  figures  in  the  chapter  comparisons  are  made 
between  theoretical  and  experimental  waveforms.   The  source  waveform  used 
in  the  theoretical  computations  was  the  observed  source  waveform.   Hence, 
the  theoretical  and  experimental  waveforms  at  the  point  x=l  ft  are 
identical. 
2.   Measurement  of  Thermoviscous  Attenuation  Coefficients 

Figure  (6-1 )  shows  the  measured  attenuation  in  the  2  in.  i.d. 
tube  as  a  function  of  frequency.   Each  measurement,  except  the  500  Hz 
measurement,  was  repeated  several  times.   The  data  is  represented  in  the 
figure  by  a  point  at  the  measurement  mean  value  and  a  bracketed  line 
representing  the  standard  deviation  of  the  measurement.   The  500  Hz 
measurement  was  only  performed  once.   Hence,  we  do  not  have  a  standard 
deviation  for  it. 

The  solid  line  in  Fig.  (6-1 )  is  a  theoretical  curve  obtained 
by  use  of  Eq.  (3-64)-   Our  measured  values  are  approximately  10$  greater 
than  the  corresponding  theoretical  values.   The  deviation  between  our 
experimental  results  and  the  Kirchhoff  theory  is  similar  to  other  reported 
deviations,  e.g.,  8$  by  Pernet  and  Payne  (1969)  and  13$  by  Mariens  (1957)  • 
The  values  used  to  evaluate  Eq.  (3-64)  were  those  for  72°F,  50$  relative 
humidity,  and  760  mm  Hg.   Handbook  values  for  these  environmental 
conditions  are  c  =  3A63I  cm/sec,  p  =  1.191  (10   )  g/cm  , 
a  =   1.85  (10"  )  g/sec-cm,  and  P  =  0.71. 
3'   Results  Concerned  with  Waves  That  are  Sinusoidal  at  the  Source 

The  basic  purpose  of  our  theoretical  and  experimental 
investigation  of  the  propagation  of  an  originally  sinusoidal  signal  was 


-     0.04 
Z     0.03  \ 


i 


u 


0.02 


128 


1 — I — I    I  I  I 


LU 

o 

u     0.01 

z 

O  0.008  f 

3  0006 
2  0.005 


<  0.004 


200 


I L 


400      600         1000  2000 

FREQUENCY  -  Hz 


J ■    i    ■ 

4000 


FIGURE  6-1 

A  COMPARISON  OF  KIRCHHOFF  ATTENUATION 

THEORY  (SOLID  LINE)  WITH  OBSERVED  DATA 


129 


to  evaluate  the  performance  of  our  computer  algorithm  and  our  equipment. 
This  section  then  is  devoted  to  the  results  of  our  developmental  testing. 
The  experimental  results  presented  here  are  similar  to  those  of 
McKittrick  et  al.  (1967)  and  some  of  those  of  Pernet  and  Payne  (1969). 
Thus  we  cannot  claim  any  novelty  as  to  the  measurements.   We  can, 
however,  claim  originality  for  our  theoretical  results.   Modified  weak- 
shock  theory  is  more  general  than  previous  theories  dealing  with  the 
propagation  of  periodic  signals  in  pipes. 

In  Figs.  (6-2),  (6-3),  (6-4),  and  (6-5)  experimentally  observed 
waveforms  are  compared  to  waveforms  computed  using  modified  weak-shock 
theory.   The  fundamental  frequency  attenuation  coefficient  used  in  the 
computation  algorithm  was  obtained  from  the  experimental  data  shown  in 
Fig.  (6.-1 ).   Basic  information  pertaining  to  the  waveforms  shown  in 
Figs.  (6-2)  through  (6-5)  is  summarized  in  Table  V. 

TABLE  V 

INFORMATION  PERTINENT  TO  THE  WAVEFORMS 
SHOWN  IN  FIGS.  (6-2)  THROUGH  (6-5) 


Figure 

Fundamental 
Frequency 
(Hz) 

Source  SPL 
(dB) 

Sh 

ock  Formation 
Distance* 
(ft) 

6-2 

800 

159.2 

10.42 

6-3 

2000 

160.0 

3.80 

6-4 

3500 

148.8 

8.26 

6-5 

2000 

160.0 

3.80 

*  The  shock  formation  distance  x  shown  in  Table  V  was  calculated  using  the 
relation 


x 


i-1 


(pek)"-1    .  (6-1) 


Equation  (6-1 )  is  derived  using  dissipationless  nonlinear  theory.   Since 
we  include  dissipation  in  our  modified  weak-shock  theory,  the  actual  shock 
formation  distances  in  our  experiments  were  somewhat  greater  than  the 
distances  shown  in  Table  V. 
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In  Fig.  (6-2)  we  present  a  comparison  of  theoretical  and 
experimental  waveforms  at  five  different  locations.   The  fundamental 
frequency  is  800  Hz.   The  dashed  waveforms  were  calculated  by  omitting 
the  wall  effects  portion  of  the  algorithm.   Thus  the  dashed  waveforms  are 
the  waveforms  that  one  would  expect  for  plane  waves  in  an  open  medium 
and,  of  course,  are  the  same  waveforms  one  could  compute  analytically 
from  ordinary  weak-shock  theory,  if  the  source  waveform  can  be  described 
analytically.   The  source  waveform  is  already  somewhat  distorted,  primar- 
ily because  of  nonlinear it ies  in  our  sound  source.   The  source  distortion 
does  not  concern  us,  however,  because  our  interest  is  in  the  change  in 
the  waveform  given  some  initial  state. 

The  observed  waveforms  agree  closely  with  those  predicted 
using  modified  weak-shock  theory  for  all  measurement  distances  reported. 
At  the  point  x=13  ft  the  observed  wave  is  still  fairly  symmetric.   The 
only  apparent  effect  of  the  tube  wall  on  the  observed  signal  is  a  small 
decrease  in  the  amplitude  over  that  predicted  using  ordinary  weak-shock 
theory.   The  effects  of  dispersion  are  not  evident  at  this  point  because 
dispersion  is  a  cumulative  effect  and  the  wave  has  not  propagated  a  great 
enough  distance  for  it  to  be  evident.   At  the  point  x=37  ft  and  beyond 
the  observed  waveforms  exhibit  the  asymmetry  previously  described. 
Notice  that  the  wave  peak  is  rounded,  the  trough  is  sharply  angled,  and 
the  shock  is  to  the  left  of  the  zero  crossing  at  the  source.   The  minor 
"raggedness"  evident  on  the  experimental  oscillograms  will  be  considered 
in  our  discussion  of  Fig.  (6-5). 

What  is  the  cause  of  the  observed  waveform  asymmetry?  In  order 
to  answer  this  question  let  us  assume  for  the  moment  that  there  is  no 
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dispersive  or  attenuative  effect  due  to  the  pipe  wall.   In  that  case  the 
small-signal  sound  speed  for  a  sinusoid  of  any  frequency  propagating  in 
the  pipe  is  the  free  field  sound  speed  c  .   Now  let  us  assume  further  that 
a  finite -amplitude  wave  is  propagating  in  the  pipe.   Nonlinear  effects 
cause  the  generation  of  harmonic  frequencies.   In  the  absence  of  dispersion 
each  harmonic*  propagates  at  the  speed  c  .   Hence  the  various  harmonics 
remain  in  step.   Eventually  the  distortion  of  the  waveform  by  nonlinear 
effects  leads  to  a  sawtooth  wave.   Analysis  of  the  symmetric  shape  and 
well  known  Fourier  series  for  the  sawtooth  amply  demonstrates  that  the 
various  harmonics  are  indeed  in  step.   Next  let  us  introduce  attenuation 
due  to  the  pipe  wall.   In  this  case  each  harmonic  is  attenuated  in 
accordance  with  Eq.  (3-64).   The  sawtooth  shape  will  be  altered  somewhat, 
but  the  resulting  waveform  will  still  be  symmetric  because  the  harmonics 
will  have  remained  in  step.   Let  us  finally  admit  the  presence  of  disper- 
sion.  Now  each  harmonic  travels  at  a  speed  given  by  Eq.  (3-63 )•   Since 
higher  frequency  harmonics  travel  faster  than  lower  frequency  harmonics 
the  harmonics  cannot  stay  in  phase.   The  loss  of  phase  coherence  of  the 
various  harmonics  leads  us  to  conclude  that  dispersion  is  the  cause  of 
the  waveform  asymmetry. 

Now  let  us  consider  the  particular  shape  of  the  asymmetric 
waveforms.   The  fact  that  higher  frequency  harmonics  travel  faster 
than  the  lower  frequency  ones  means  that  the  troughs  of  the  various 


We  draw  an  important  distinction  here  between  the  propagation  velocity 
of  a  wavelet  and  the  propagation  velocity  of  some  harmonic.   The  wavelet 
travels  at  a  speed  given  by  use  of  Eq.  (2-4).   The  harmonics  are  being 
continuously  generated  by  nonlinear  effects  but  the  propagation  velocity 
of  the  harmonic  itself  is  given  in  the  absence  of  dispersion  by  the  small- 
signal  sound  speed. 
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harmonics  will  become  lined  up  whereas  the  peaks  will  spread  apart. 
Therefore  the  trough  of  the  overall  waveform  will  be  sharply  angled 
whereas  the  peak  will  be  broad  and  rounded;  at  the  same  time  the  zero 
crossing  will  shift  to  the  left  of  its  original  position.   The  reader 
may  quickly  convince  himself  that  our  analysis  is  plausible  by  sketching 
the  first  several  harmonics  of  a  sawtooth  and  then  displacing  the  various 
harmonics  for  the  case  in  which  speed  increases  with  frequency.   The 
resulting  overall  waveform  will  exhibit  the  asymmetry  described  here. 

Figure  (6-3 )*  shows  a  comparison  of  2000  Hz  fundamental  frequency 
experimental  and  theoretical  waves.   The  dashed  weak-shock  waveforms 
included  in  Fig.  (6-2)  have  been  omitted  in  Fig.  (6-3).   At  the  point 
x=13  ft  the  wave  has  traveled  about  3  shock  formation  lengths  so  that 
on  the  basis  of  ordinary  weak-shock  theory  we  would  expect  a  sawtooth 
wave.   Because  of  wall  effects,  however,  there  is  some  rounding  of  the 
peak  of  the  wave. 

In  Fig.  (6-4)  the  fundamental  frequency  of  the  waves  is  3500  Hz. 
The  purpose  of  Fig. (6- 4)  is  to  demonstrate  that  our  method  is  applicable 
to  waves  having  fundamentals  at  least  as  high  as  3500  Hz.   We  have  not 
quantitatively  defined  the  point  at  which  a  finite-amplitude  wave  no  longer 
contains  shocks. *~*  Qualitatively,  however,  we  can  observe  that  the  shocks 
in  the  waveforms  at  the  point  x=85  ft  are  badly  degraded. 


Figure  (o-3)  has  previously  been  reported  by  Pestorius  and  Blackstock 
(1972). 

Blackstock  (1966)  calculates  the  inception  of  "old  age"  as  occurring 
at  the  point 

x  =  aQ-  (6-2) 

In  the  case  shown  in  Fig.  (6-4)  evaluation  of  Eq.  (6-2)  yields  a  ~*-=kO   ft. 
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In  Fig.  (6-5)  the  sets  of  waveforms  at  the  points  x=13  ft, 
x=3T  ft,  x=6l  ft,  and  x=85  ft  shown  in  Fig.  (6-3)  have  been  enlarged  and 
overlaid. .  The  agreement  between  theory  and  experiment  in  each  compari- 
son is  excellent.   In  particular  note  that  the  waveforms  agree  at  the 
point  x=85  ft  as  well  as  they  agree  at  the  shorter  distances,  thereby 
indicating  that  there  are  no  apparent  cumulative  errors  in  our  method. 
Similar  comparisons  of  the  waveforms  shown  in  Figs.  (6-3)  and  (6-h) 
have  been  made.   In  general  the  experimental  and  theoretical  waveforms 
shown  in  Figs.  (6-2)  and  (6-3)  agree  with  one  another  in  a  manner  similar 
to  the  agreement  shown  in  Fig.  (6-5). 

The  minor  "raggedness"  in  the  experimental  data  appears  to  be 

due  to  scattering  from  minor  irregularities  in  the  pipe  and  pipe  flanges 

near  the  microphone.   The  raggedness  appears  to  be  shock  excited  since  it 

is  not  apparent  until  shocks  are  present.   The  pattern  of  these  waveform 

irregularities  is  very  repeatable  cycle  to  cycle,  but  the  pattern  can 

be  varied  at  any  particular  measurement  location  by  disassembling  and 

reassembling  the  pipe  flange  at  that  location.   We  have  not  conducted 

any  further  analysis  of  the  waveform  irregularities  since,  no  matter 

what  their  origin  is,  they  are  minor  compared  to  the  effects  that  we  are 

studying. 

h.      A  Comparison  of  Results  Obtained  Using  Waves  That  are  Sinusoidal 
at  the  Source  with  Those  of  Other  Researchers 

In  Figs.  (6-6)  and  (6-7 )  we  compare  some  of  our  results 

obtained  using  waves  that  are  originally  sinusoidal  with  those  of  Pernet 

and  Payne  (1971a)  and  Coppens  (igjl) ,   respectively.   In  Fig.  (6-8)  we 

compare  results  obtained  using  weak-shock  theory,  modified  weak-shock 


138 


theory,  and  the  theory  derived  by  Rudnick  (1952).   Each  of  these 
comparisons  will  be  discussed  in  turn. 

Pernet  and  Payne  (1971 )  were  concerned  with  wave  propagation 
problems  in  which  finite -amplitude  distortion  is  significant  but  shocks 
do  not  form.   In  the  developmental  testing  of  their  analytic  model, 
high  intensity  (e.g.,  160  dB)  sinusoids  were  generated  and  the  levels 
of  the  fundamental  and  the  first  nine  harmonics  recorded  as  a  function 
of  distance.   Although  the  intensity  of  these  sinusoids  was  too  great 
for  analysis  using  their  method,  they  did  given  the  experimental  data  for 
a  500  Hz  sinusoid  and  its  first  nine  harmonics.   We  have  modified  program 

ATENDIS5*  in  such  a  way  that  we  may  calculate  the  amplitude  of  the 
fundamental  and  the  2nd  through  5th  harmonics  of  a  500  Hz  sinusoid  as 

a  function  of  distance.   In  Fig.  (6-6)  we  have  compared  the  Pernet  and 
Payne  data  (solid  curves  )  for  these  signals  with  our  theoretical  predic- 
tions (symbols).   In  the  calculation  of  the  theoretical  data,  our  measured 
value  of  the  fundamental  frequency  attenuation  coefficient  was  used 
[see  Fig.  (6-1)].  The  agreement  between  the  experimental  data  and  theory 
is  quite  good  although  our  theory  points  lie  slightly  above  the  data. 
Of  course,  since  we  do  not  know  the  accuracy  of  their  measurements  an 
absolutely  quantitative  comparison  is  difficult  to  make.   The  important 
fact  is  that  the  Pernet  and  Payne  data  provide  an  independent  check  of 
modified  weak-shock  theory  as  implemented  by  our  algorithm.   The  Pernet 


•x- 
For  those  readers  who  did  not  read  Chapter  4,  program  ATENDIS5  is 

a  computer  program  wherein  we  combine  weak-shock  theory  and  tube  wall 

attenuation  and  dispersion.   Program  ATENDIS5  was  used  to  compute  the 

theoretical  results  presented  in  Sections  3  and  4.   A  similar  program, 

Program  NSATDIS,  was  used  to  compute  the  results  in  Sections  5,  6, 

and  7.   Programs  ATENDIS5  and  NSATDIS  are  very  similar  in  content. 

The  reason  for  using  two  separate  programs  was  primarily  for  convenience, 
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and  Payne  paper  contains  several  other  sets  of  data  that  could  be 
checked  against  our  theory.   These  checks  have  not  been  performed  for 
want  of  time. 

Coppens'  (l9Tl)  results  are  purely  theoretical  and  are  valid 
only  for  propagation  distances  that  preclude  shock  formation.   Several 
plots  are  given  that  depict  the  waveform  distortion  that  occurs  as  an 
originally  sinusoidal  wave  propagates  in  a  pipe.   These  plots  are  similar 
in  content  to  the  theoretical  curves  shown  in  Figs.  (6-2)  through  (6-5). 
Coppens  defined  two  dimensionless  parameters  as  follows: 

I   =  ax   ,  (6-3) 

and 

Mp   pek  (c   ,  x 

1     0 

Evaluation  of  Eqs.  (6-3)  and  (6-k)   at  the  point  x=6l  ft  using  a  source 

sinusoid  having  a  sound  pressure  level  value  of  1^2  dB  and  a  frequency  of 

2000  Hz  yields  the  values  ,2=0.  99  and  Mp/S  =1.0  .   These  two  values  closely 

match  the  parameter  values  for  one  of  Coppens'  published  waveforms. 

Hence,  we  may  directly  compare  our  theory  and  experimental  data  with 

Coppens1  theory.   In  Fig.  (6-7 )  we  have  compared  modified  weak-shock  theory 

(solid  line),  our  experimental  data  points  (circles ),  and  Coppens'  theory 

(dotted  line).   The  agreement  between  our  theory  and  data  is  excellent. 

Coppens'  theory  also  agrees  closely  with  ours.   The  differences  between 

the  two  theoretical  curves  are  probably  due  to  the  difference  in  the  value 

of  the  fundamental  frequency  attenuation  coefficient  used  in  each  analysis. 
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Coppens  used  the  Kirchhoff  theory  value  while  we  used  a  value  that  we 
had  measured.   The  measured  value  exceeded  the  theoretical  value  by 
10$  and  the  difference  is  probably  the  cause  of  the  greater  signal 
attenuation  evident  in  our  results. 

In  modified  weak-shock  theory  we  combine  shock  effects  and 
pipe  wall  effects  in  an  iterative  fashion.   A  possible  alternate  compu- 
tational scheme  consists  of  combining  the  nonlinear  and  pipe  wall  attenu- 
ation (but  not  dispersion)  in  a  nonlinear  differential  equation.   Rudnick 
(1952)  derived  such  an  equation  for  the  propagation  of  sawtooth  waves. 
Recasting  Rudnick' s  equation  in  terms  of  the  variables  used  here,  we  obtain 

dP 

_^£=  _aP    __gk     P   2 

dx       max        2   max  ( 6-5 ) 

up   c  v  ■" 

0  0 

where 

P    =  peak  acoustic  pressure,  and 

max  ' 

a   =  the  attenuation  coefficient  calculated  using  Rudnick' s 
"zero-order  correction"  method. 
The  first  term  on  the  right  hand  side  of  Eq.  (6-5)  accounts  for  thermo- 
viscous  attenuation  and  the  last  term  accounts  for  attenuation  due  to  the 
effects  of  nonlinear  propagation.   The  solution  of  Eq.  (6-5)  subjects  the 
boundary  condition  that  P    =  ?        (0)   at  x=0  is 


Pmax(0)  e-ax 

P       =  TTTfr 7-r-T .  (6-6) 

max       pkP   (0)        ^ 
..  |    max    /,    -Ox  x 
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In  Fig.  (6-8)  we  compare  our  theory  and  experimental  results 
with  predictions  obtained  using  Eq.  (6-6).   Referring  to  the  waveforms 
presented  in  Fig.  (6-3),  we  selected  the  experimental  waveform  at  the 
point  x=13  ft  as  a  starting  point  for  our  calculations.   This  waveform 
is  appropriate  because  it  is  almost  sawtoothed  in  shape  (use  of  Eq.  (6-6) 
requires  that  the  wave  be  a  sawtooth  to  begin  with).   The  data  presented 
in  Fig.  (6-8)  shows  the  change  in  peak-to-peak  amplitude  of  the  wave 
as  it  propagates.   Curve  1  is  the  result  obtained  using  weak-shock  theory, 
curve  2  is  the  result  obtained  using  modified  weak-shock  theory  and 
curve  3  is  "the  result  of  calculations  based  on  Eq.  (6-6).   The  experimental 
points  confirm  the  middle  curve.   Use  of  weak-shock  theory  predicts  too 
little  attenuation  because  absorption  due  to  the  tube  wall  boundary  layer 
is  neglected.   Yet  when  we  linearly  add  the  wall  absorption  we  predict 
too  much  attenuation.   The  key  to  the  agreement  between  our  theory  and 
experimental  data  is  our  inclusion  of  dispersion  in  modified  weak-shock 
theory.   The  severe  rounding  of  the  peaks  of  the  shocks  caused  by  disper- 
sion in  effect  reduces  the  shock  amplitude.   Referring  to  Eq.  (3-17), 
we  see  that  attenuation  at  the  shock  depends  upon  the  shock  amplitude. 
Hence,  dispersion  inhibits  the  attenuation  associated  with  the  shock. 
The  combined  effects  of  the  reduced  shock  attenuation  and  the  tube  wall 
attenuation  then  yield  the  results  shown  in  the  middle  curve. 


Figure   (6-8)  has  previously  been  reported  by  Pestorius  and  Blackstock 
(1972). 
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5.   Noise  Waveform  Results 

As  we  noted  in  the  introductory  chapter,  we  are  especially 
interested  in  the  finite-amplitude  distortion  of  noise  waveforms.   In 
Figs.  (6-9)  and  (6-11 )  comparisons  are  made  between  measured  and  computed 
noise  waveforms.   In  Fig.  (6-10 )  comparison  is  made  between  weak-shock 
theory  and  modified  weak-shock  theory  using  some  of  the  waveforms  in 
Fig.  (6-9).   The  theoretical  results  in  Figs.  (6-9),  (6-10 ),  and  (6-11 ) 
were  computed  using  program  NSATDIS.   As  was  done  in  the  examples  in 
Section  3,  the  source  waveform  for  use  in  the  noise  waveform  computations 
was  obtained  by  sampling  the  observed  wave  at  the  point  x=l  ft.   In  these 
computations  the  theoretical  value  of  Oi     was  used  rather  than  the  measured 
value  because  a  direct  measurement  could  not  be  made  at  the  fundamental 
frequency  [6l.h   Hz  for  the  waves  shown  in  Fig.  (6-9)  and  78.2  Hz  for  the 
waves  shown  in  Fig.  (6-ll)]. 

Computed  and  measured  waveforms  are  shown  in  Fig.  (6-9).   The 
SPL  at  the  source  was  160  dB.   The  agreement  between  theory  and  experiment 
is  quite  good.   Note,  however,  that  unlike  the  results  given  in  Section  3, 
the  apparent  agreement  between  the  theoretical  results  shown  and  the 
measurements  deteriorates  slightly  as  the  propagation  distance  gets  large. 
We  believe  that  adjustment  of  the  incremental  distance  steps  in  our 
program  could  minimize  this  small  disagreement.   However,  the  necessary 
adjustment  was  not  deemed  worthwhile  in  terms  of  the  cost  of  the  additional 
computer  time  required  and,  hence,  was  not  performed. 


Figures  (6-9)  and  (6-10 )  have  previously  been  reported  by  Pestorius  and 
Blackstock  (1975a). 
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In  Fig.  (6-10)  the  computed  -waveforms  at  1J5  ft  and  85  ft,  shown 
in  Fig.  (6-9) ,    are  repeated  and  compared  to  waveforms  computed  by  using 
the  algorithm  based  on  weak-shock  theory  alone.   It  is  clear  that  the  wall 
effects  have  a  very  significant  influence  on  the  waveform,  particularly  at 
the  greater  distance.   This  increased  influence  is  caused  by  the  fact  that 
the  effects  of  attenuation  and  dispersion  are  cumulative  with  distance 
and  the  fact  that  the  nonlinear  effects  grow  weaker  with  distance. 

Figure  (6-11 )  is  similar  to  Fig.  (6-9)  except  that  a  different 
sample  pulse  is  shown.   Figure  (6-11 )  is  included  to  demonstrate  that  the 
results  shown  in  Fig.  (6-9)  are  not  atypical.   The  sound  pressure  level 
for  the  pulse  shown  in  Fig.  (6-11 )  was  155  dB. 

The  waveforms  shown  in  Figs.  (6-9)  and  (6-11 )  undergo  marked 
change  as  they  propagate.   The  source  waveforms  in  both  examples  are 
free  of  shocks.   At  the  point  x=13  ft,  however,  the  waveforms  in  both 
examples  have  numerous  shocks.   Since  individual  positive  sloping 
sections  are  very  steep,  the  first  shocks  form  within  a  few  feet  of  the 
source  point. 

Let  us  now  restrict  our  attention  to  the  results  shown  in 
Fig.  (6-9).   At  the  points  x=13  ft  and  x=3T  ft  "the  waves  are  still  very 
intense  and  the  rounding  effect  at  the  wave  peaks  due  to  dispersion  is 

not  very  evident.   The  lack  of  noticeable  rounding  indicates  that  the 
dispersive  effects  are  being  obscured  by  the  much  stronger  nonlinear 
effects.   At  the  points  x=6l  ft  and  x=85  ft  the  rounding  of  the  wave 
peaks  is  very  evident.   Note  also  that  the  propagating  wave  simplifies 
considerably  as  the  shocks  merge  with  one  another.   These  waveform 
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changes  can,  of  course,  be  used  to  predict  changes  in  the  wave  spectrum 
as  a  function  of  distance.   We  shall,  however,  defer  the  discussion  of 
spectra  until  Section  6. 

In  the  preceding  paragraph  we  restricted  our  attention  to  the 
results  shown  in  Fig.  (6-9),  but,  in  general,  the  above  remarks  apply  also 
to  the  results  shown  in  Fig.  (6-11 ).   The  source  SPL  for  this  latter 
example  was  lower  (155  dB).   Therefore  the  effects  of  dispersion  become 
evident  sooner  than  in  the  example  shown  in  Fig.  (6-9) • 
6.   Predictions  and  Measurements  of  Noise  Spectra 

In  this  section  we  present  three  sets  of  spectral  results. 
The  first  two  sets  are  computer  generated  and  include  the  effects  of 
nonlinear  propagation,  but  do  not  include  the  effects  of  the  tube  wall. 
The  last  group  of  spectra  were  obtained  experimentally. 

We  used  the  waveforms  at  the  point  x=l  ft  in  Figs.  (6-9)  and 
(6-11 )  as  the  input  data  for  program  DCTTALL3*  in  order  to  compute  the 
first  set  of  results.   Spectra  were  computed  at  the  point  x=l  ft,  and, 
following  "computer  propagation",  at  the  points  x=U9  ft  and  x=73  ft* 
These  results  appear  in  Figs.  (6-12 )*•*  and  (6-13 ),  respectively. 

In  the  second  group  of  results  the  spectra  are  entirely  computer 
generated.   The  source  waveforms  were  obtained  using  the  digitally  filter- 
ed output  of  a  computer  random  number  generator.   The  results  of  these 
computations  appear  in  Fig.  (6-lU).   As  in  the  previous  group  of  results 


For  those  who  did  not  read  Chapter  h,   D0ITALL3  is  a  computer  program 
used  to  calculate  spectra.   The  computation  method  used  is  the  Blackman 
and  Tukey  (1959)  indirect  method. 

■3Ht 

Figures  (6-12),  (6-lU),  and  (6-15 )  have  been  previously  reported  by 
Pestorius  and  Blackstock  (1973b). 
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spectra  at  the  points  x=l  ft,  x=^9  ft,  and  x=73  ft  are  presented. 

The  third  set  of  results  are  shown  in  Fig.  (6-I5).   These  results 
were  obtained  using  the  General  Radio  Wave  Analyzer  and  Recorder.   A 
50  Hz  crystal  was  used  in  the  wave  analyzer.   The  frequency  range  was 
scanned  very  slowly  to  give  a  long  integration  time.   Stable  statistical 
estimates  were  obtained  in  this  manner.   The  noise  signal  was  produced  by 
bandpass  filtering  the  output  of  the  ELGENCO  Model  603A  noise  generator, 
power  amplifying  this  output,  and  applying  it  to  a  horn  driver  attached 
to  the  plane  wave  pipe. 

We  will  now  consider  in  turn  each  of  these  sets  of  results.   In 
Figs.  (6-12)  and  (6-I3)  three  spectra  are  presented.   The  equivalent  band- 
width of  the  spectral  estimates  is  50  Hz.   The  spectra  shown  in  both 
figures  were  normalized  by  dividing  each  spectrum  data  point  by  the  maxi- 
mum value  in  the  respective  source  spectrum.   Had  the  propagation  been 
in  accordance  with  linear  acoustic  theory,  the  three  spectra  in  each 
figure  would  have  been  identical.   The  flattening  and  spreading  of  the 
spectra  shown  in  Figs.  (6-12)  and  (6-I3)  are  the  result  of  nonlinear 
effects.   In  particular,  the  spectra  at  the  points  x=l+9  ft  and  x=73  ft 
actually  cross  over  the  source  spectrum  indicating  a  net  increase  in 
energy  at  frequencies  above  the  crossover  frequency.   In  order  to  investi- 
gate the  effects  of  aliasing  inherent  in  spectral  computations,  the 
frequency  at  which  the  source  waveforms  were  sampled  was  increased  to 
196  kHz  and  the  spectral  computations  were  repeated.   The  analysis  showed 
that  aliasing  was  minimal  and,  in  particular,  that  no  large  line  components 
were  present  above  the  frequencies  shown  in  Figs.  (6-12)  and  (6-I3). 
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The  purpose  of  this  group  of  results  then  was  to  l)  predict  the  spectral 
changes  due  to  nonlinear  effects  for  particular  noise  waveforms,  and 
2)  to  investigate  the  effects  of  aliasing  present  in  our  spectral  compu- 
tations. 

The  purpose  of  the  spectral  measurements  shown  in  Fig.  (6-lU) 
was  to  predict  as  closely  as  possible  the  actual  average  spectral  shape 
of  noise  waves  propagated  in  our  pipe.   The  3  dB  down  points  used  in  the 
digital  filter  program  were  TOO  Hz  and  2200  Hz  and  the  slopes  of  the 
skirts  of  the  source  spectrum  closely  match  those  observed  at  the  horn 
driver  output  during  our  initial  system  calibration.   The  spectra  in 
Fig.  (6-14)  are  smoother  than  those  in  Figs.  (6-12)  and  (6-13)  because  of 
increased  averaging. 

In  Fig.  (6-I5)  actual  experimental  spectra  are  shown.   The 
results  are  very  similar  to  the  results  shown  in  Figs.  (6-12),  (6-13), 
and  (6-lU).   In  particular,  the  high  frequency  crossover  shown  in  all 
three  of  the  former  figures  appears  in  Fig.  (6-I5).   The  low  frequency 
crossover  predicted  in  the  results  shown  in  Fig.  (6-14)  was  observed  and 
appears  in  Fig.  (6-I5).   Since  the  source  spectrum  shown  in  Fig.  (6-lU) 
was  tailored  to  match  our  average  source  conditions  as  closely  as  possible, 
the  verification  of  the  low  frequency  crossover  prediction  is  especially 
significant.   The  total  dissipation  suffered  by  spectral  components  shown 
in  Fig.  (6-I5)  is  greater  than  that  suffered  by  the  same  components  shown 
in  the  computed  results  because  the  spectra  shown  in  Fig.  (6-15/  have 
suffered  tube  wall  losses  as  well  as  nonlinear  attenuation.   The  local 
spectral  peaks  at  about  8  kHz  shown  in  Fig.  (6-I5)  are  apparently  due  to 
a  resonance  In  the  horn  driver. 
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Let  us  now  consider  the  change  in  the  wave  spectra  interpreted  in  terms 
of  the  waveform  changes  noted  in  the  previous  section.  We  observed  that 
the  waves  simplified  significantly  as  a  function  of  distance.   In  particu- 
lar the  number  of  zero  crossings  diminishes.   The  waveform  simplification 
represents  a  transfer  of  energy  to  the  lower  end  of  the  spectrum. 
Conversely,  the  presence  of  nonlinear  steepening  means  that  the  high 
frequency  end  of  the  spectrum  is  also  enhanced.   In  fact,  the  high  fre- 
quency build-up  commences  as  soon  as  the  wave  is  formed  whereas  the  low 
frequency  effect  requires  the  presence  of  shock  fronts.   Once  the  shocks 
in  a  particular  noise  sample  have  all  formed,  the  flow  of  energy  will 
proceed  only  toward  the  low  end  of  the  spectrum.   Thus  the  overall  effect 
of  nonlinear  propagation  is  to  flatten  out  the  wave  spectrum. 

Alternatively,  these  same  effects  can  be  understood  as  the  result 
of  the  intermodulation  of  various  source  spectral  components.   If  we  think 
of  the  source  spectrum  as  a  band  of  closely  spaced  line  components,  each 
component  will  interact  with  each  other  component  to  produce  difference 
frequencies  (low  frequency  build-up)  and  sum  and  harmonic  frequencies 
(high  frequency  build-up).   These  two  explanations,  one  in  terms  of  time 
domain  behavior  and  the  other  in  terms  of  the  frequency  domain,  are,  of 
course,  equivalent. 

The  noise  results  presented  here  and  in  Section  5  represent  the 
heart  of  the  results  presented  in  this  work.   In  the  next  few  paragraphs 
we  will  note  the  most  significant  aspects  of  these  results. 

Let  us  consider  the  waveform  results  first.   Nonlinear  effects 
cause  a  systematic  simplification  of  noise  waveforms  as  a  function  of 
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distance.   In  the  two  examples  shown  in  Section  5  the  waves  approach  a 
shape  consisting  of  shock  fronts  separated  by  negatively  sloping  sections. 
Because  of  dispersion  these  sections  are  somewhat  curved.   If  we  draw 
straight  line  tangents  to  these  sections  we  find  that  all  the  tangent 
slopes  are  very  similar.   Hence  (ignoring  the  effects  of  dispersion)  we 
observe  experimentally  the  asymptotic  waveform  predicted  by  Burgers. 

Let  us  now  turn  to  the  spectral  results.   It  is  quite  clear 
that  the  effects  of  nonlinear  propagation  cause  marked  changes  in  noise 
spectra.   In  particular  the  high  and  low  frequency  spectral  components 
may  exhibit  a  net  increase  in  energy  as  a  function  of  distance.   The 
high  frequency  build-up  commences  as  soon  as  the  finite -amplitude  wave 
is  formed.   The  low  frequency  effect  requires  the  formation  of  shocks. 
The  cascading  of  energy  to  lower  frequencies  observed  here  verifies 
Kraichnan's  (1968)  predictions  concerning  spectral  shifting  to  lower 
wave numbers. 

Since  the  low  frequency  energy  cascade  requires  the  presence 
of  shocks,  we  expect  to  observe  it  at  a  higher  SPL  than  that  required  to 
observe  the  high  frequency  cascade.   We  have  begun  a  preliminary  study 
of  the  threshold  SPL  necessary  to  cause  the  high  frequency  crossover  to 
occur.   In  our  plane  wave  tube  the  crossover  occurs  for  source  SPLs  of 
132  to  135  dB  and  above.   We  have  not,  to  date,  investigated  the  threshold 
of  the  low  frequency  crossover  except  to  note  that  it  is  not  evident  at 
135  dB. 


159 


7.   The  Relation  of  Our  Results  Obtained  Using  Noise  Waveforms  to 
Results  Obtained  Using  the  "Burgers   Turbulence"  Model 

In  Chapter  3  we  derived,  via  weak-shock  theory,  the  result 
that  for  large  value  of  x  the  shape  of  finite-amplitude  noise  waves  is 
equivalent  to  the  random  "sawtooth  profile"  obtained  using  the  Burgers 
turbulence  model.   In  particular  we  found  that  the  slopes  of  the  contin- 
uous sections  of  the  wave  became  independent  of  their  initial  slope  and 
proportional  to  x~  .   In  Fig.  (6-10)  the  waveform  at  the  point  x=85  ft 
computed  using  weak-shock  theory  clearly  exhibits  the  characteristics  of 
a  random  sawtooth.   In  addition,  if  we  ignore  the  effects  of  dispersion 
in  the  experimental  and  theoretical  results  shown  in  Figs.  (6-9)  and 
(6-11 ),  we  observe  that  the  waveforms  for  the  larger  values  of  x  look 
like  random  sawtooth  waves.   Hence,  we  have  theoretically,  via  weak-shock 
theory,  and  experimentally  confirmed  the  random  sawtooth  profile  predicted 
by  Burgers  (1972). 

In  the  spectral  results  presented  in  the  previous  section  we 
observed  both  experimentally  and  theoretically  the  waveform  simplification 
and  downward  energy  cascade  predicted  by  Kraichnan  (1968). 

We  shall  report  here  on  our  attempt  to  computationally  verify 
the  value  of  £  predicted  by  Tatsumi  and  Kida  (1972).  They  derived  the 
result  that  the  number  of  shocks  in  a  "random  sawtooth"  varies  as 


^  =  \(r-Y    ■  <2-8> 


Converting  Eq.  (2-8)  to  our  boundary  value  problem  notation,  we  obtain 
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N(x)  =  N   —      .  (6-7) 
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Tatsumi  and  Kida  predicted  for  very  large  Reynolds  number  (very  large 

amplitude  waves)  and  for  very  large  values  of  propagation  distance  that 

the  value  of  (•  should  approach  l/2.   Using  our  program  NSATDIS,  we  have 

made  a  check  of  the  validity  of  this  prediction.   Inspection  of  Figs.  (6-9) 

and  (6-ll)  show  that  at  13.  ft  the  waveforms  are  approaching  the  asymptotic 

random  sawtooth  form.   Hence,  we  selected  the  computed  waveforms  at  the 

point  x=13  ft  as  our  source  waveforms  and  13  ft  as  the  value  of  x  .   We 

0 

then  calculated,  using  the  weak-shock  theory  portion  of  program  NSATDIS, 
the  waveforms  for  values  of  x  greater  than  x  .   In  Burgers '  equation 
the  effects  of  nonlinear  steepening,  including  attenuation  at  the  shocks, 
are  combined  with  the  effects  of  thermoviscous  attenuation  in  the  contin- 
uous portion  of  the  wave.   Tatsumi  and  Kida's  use  of  the  large  Reynolds 
number  assumption,  however,  removes  the  effects  of  attenuation  in  the 
continuous  portion  of  the  wave  leaving  only  the  attenuation  at  the  shocks. 
Therefore,  our  use  of  only  the  weak-shock  theory  portion  of  program 
NSATDIS  is  consistent  with  Tatsumi  and  Kida's  analysis.  Calculations  of  £ 
were  made  using  the  waveforms  described  above.   The  source  levels  in  these 
cases,  however,  proved  to  be  too  low  to  adequately  approximate  the  assump- 
tion of  very  high  Reynolds  number.   Therefore,  using  the  computer  we 
scaled  the  waveform  shown  in  Fig.  (6-9)  to  a  source  level  of  175  dB.   The 
value  of  i   was  calculated  by  plotting  the  waveform  at  various  propagation 
distances,  counting  the  number  of  shocks  and  using  Eq.  (6-7 )•   The  results 
of  these  computations  are  presented  in  Table  VI. 
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TABLE  VI 

Calculation  of  |  for  x  =13  ft 

o 

Source  SPL  (dB)       £_ 

155  0.097 

160         0.1999 
175         0.281 


The  values  of  £  given  in  Table  VI  are  far  lower  than  the 
value  l/2  predicted  by  Tatsumi  and  Kida.   Note  that  i   appears  to  be 
increasing  as  the  source  level  is  increased  and  for  a  strong  enough  source 
might  approach  the  value  l/2.   However,  the  source  level  required  to  obtain 
the  value  1/2  would  certainly  exceed  the  levels  shown  in  Table  VI  and 
would,  therefore,  be  too  high  for  the  use  of  weak-shock  theory.   Hence, 
we  conclude  that  Tatsumi  and  Kida's  asymptotic  value  of  £  cannot  be  compu- 
tationally verified  by  the  use  of  weak-shock  theory  because  weak-shock 
theory  is  not  valid  at  the  amplitudes  required  for  £  to  approach  l/2. 


CHAPTER  7  CONCLUDING  REMARKS 

1.  Introduction 

In  Chapter  7  we  present  a  brief  summary  of  the  major  results 
of  our  research  and  draw  conclusions  based  upon  these  results.   We  then 
consider  the  general  question  of  the  importance  of  dispersion  in  the  study 
of  propagation  in  pipes.   We  conclude  by  indicating  some  possible  avenues 
of  further  experimental  and  theoretical  research. 

2.  Summary  of  Results 

The  major  results  presented  in  this  work  can  be  divided  into 
three  categories:   the  development  of  a  computational  method  for  handling 
the  propagation  in  a  tube  of  finite -amplitude  waves  of  arbitrary  source 
waveform,  results  obtained  using  originally  sinusoidal  waves,  and  results 
obtained  using  noise  waves.   We  will  summarize  in  turn  the  results  in 
each  of  these  categories. 

In  Chapters  3  and  h   we  developed  our  computational  algorithm. 
Our  method  is  based  upon  weak-shock  theory,  with  the  added  provision  that 
we  account  for  absorption  and  dispersion  due  to  the  tube  wall  boundary 
layer.   This  provision  is  very  important  because  it  enables  us  to  accurate- 
ly predict  the  shape  of  sound  waves  propagating  in  a  pipe.   In  addition, 
it  extends  the  validity  of  weak-shock  theory  downward  into  the  region  of 
waves  too  weak  to  be  treated  by  ordinary  weak  shock  theory  or  even  pertu- 
bation  analysis  but  still  too  strong  for  linear  theory  to  be  adequate. 
As  Blackstock  (1972)  notes,  ordinary  weak-shock  theory  does  not,  in 
general,  work  for  very  weak  waves  because  it  does  not  contain  any  attenu- 
ation effects  other  than  that  due  to  the  shocks.   The  fact  that  the 

modified  weak-shock  algorithm  includes  tube  wall  attenuation  as  well  as 
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shock  attenuation  means  that  our  method  should  not  have  a  lower  strength 
limit. 

Although  we  developed  the  algorithm  for  the  specific  case  of 
propagation  in  tubes,  the  extension  of  the  method  to  describe  the  propaga- 
tion of  waves  in  an  unbounded  medium  is  straightforward. 

The  studies  of  waves  generated  by  a  sinus oidally  driven  source 
were  carried  out  to  test  our  theoretical- computational  algorithm.   In 
Figs.  (6-2)  through  (6-5)  we  see  excellent  agreement  between  our  theoretical 
predictions  and  observed  data.   Comparison  with  the  measurements  made  by 
Pernet  and  Payne  (1971a)  provides  an  independent  verification  of  our 
method.   The  comparison  with  Coppens'  (1971 )  theoretical  results  gives 
assurance  that  our  method  is  valid  for  waves  too  weak  to  form  shocks 
as  well  as  for  stronger  ones.   We  cannot  claim  novelty  for  the  experimen- 
tal results  obtained  using  sinusoidal  waves  because  the  waveform  distortion, 
including  the  asymmetry,  has  been  reported  previously  by  McKittrick  et  al. 
(1967)  and  Pernet  and  Payne  (1969). 

As  we  stated  in  Chapter  1,  our  major  interest  has  been  in  the 
propagation  of  finite-amplitude  noise.  In  Chapter  6  two  types  of  noise 
results  are  presented:   waveforms  and  spectra.   The  former  show  very 
clearly  the  changes  that  take  place  in  the  waveform  because  of  nonlinear ity, 
for  example,  waveform  steepening,  reduction  in  the  number  of  zero  crossings, 
and  general  waveform  simplification.   Two  examples  of  noise  pulses  are 
presented  in  Chapter  6  in  Figs.  (6-9)  and  (6-ll).   In  both  cases  agreement 
between  theory  and  experiment  is,  in  general,  excellent.   Analysis  of  the 
results  shown  in  Fig.  (6-10 )  indicate  the  importance  of  the  inclusion  of 
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the  tube  wall  attenuation  and  dispersion  in  our  theory.   The  results  shown 
in  Figs.  (6-9),    (6-10),  and  (6-11 )  amply  illustrate  the  effects  of  non- 
linear propagation. 

The  changes  in  the  noise  waveforms  with  distance  are  reflected 
in  the  waveform  spectra  shown  in  Chapter  6.   Examples  are  presented  in 
Figs.  (6-12),  (6-13),  (6-14),  and  (6-15 ).   In  Figs.  (6-12)  and  (6-13) 
computed  spectra,  calculated  using  the  source  waveforms  in  Figs.  (6-9) 
and  (6-11 ),  respectively,  are  shown.   The  spectra  in  Fig.  (6-14)  are 
entirely  computer  generated.   In  Fig. (6-I5)  experimental  spectra  obtained 
using  a  wave  analyzer  are  depicted.   Comparison  of  Figs.  (6-14)  and  (6-I5) 
demonstrates  general  agreement  between  theory  and  experiment.   The  spectra 
in  Figs.  (6-12 )  through  (6-I5)  are  perhaps  our  most  important  results. 
Analysis  of  these  figures  shows  that  the  spectra  associated  with  both 
experimental  and  computed  waves  undergo  marked  change  as  the  wave  propa- 
gates.  These  changes  have  not  been  observed  by  previous  researchers 
although  Kraichnan  (1968)  and  Ffowcs  Williams  (1965)  did  predict  the  low 
frequency  energy  cascade.   The  high  and  low  frequency  cascades  and  actual 
spectral  crossovers  mean  that  nonlinear  effects  may  cause  a  net  increase 
in  energy  at  some  frequencies  as  the  waves  propagate. 

In  Chapter  6,  Section  J,   we  reviewed  results  that  showed 
similarities  between  results  obtained  using  weak-shock  theory  and  those 
obtained  using  Burgers'  equation.   In  addition  we  attempted  to  verify 
Tatsumi  and  Kida's  (1972)  asymptotic  value  of  the  quantity  |.   The  value 
of  £  is  a  measure  of  the  randomness  of  a  turbulent  field,  or,  as  in 
our  case,  a  noise  wave.   We  were  unable  to  verify  the  predicted  value. 
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3.   Conclusions 

We  will  consider  first  our  theoretical-computational  method. 
We  conclude  that  our  modified  weak-shock  algorithm  is  novel  in  that  it  is 
more  general  than  previously  proposed  schemes  that  deal  with  the  propa- 
gation of  plane  finite-amplitude  waves  in  pipes.   It  is  original  in  that 
results  produced  using  it  are  valid  both  before  and  after  shock  formation, 
the  results  include  both  nonlinear  and  pipe  wall  effects,  and  the  method 
works  for  arbitrarily  shaped  input  waves.   Previous  methods,  while 
containing  some  of  the  above  features,  do  not  contain  them  all. 

We  will  consider  next  results  obtained  using  a  sinusoidally 
driven  source.   Our  theoretical- computational  model  has  been  experimen- 
tally verified  when  using  waves  that  are  sinusoidal  at  the  source  over 
the  frequency  range  of  500  to  3500  Hz  and  the  source  sound  pressure  level 
range  of  142  to  160  dB.   Although  the  actual  experimental  verification 
of  the  theory  is  at  present  limited  to  the  frequency  and  source  level 
ranges  reported  here,  we  believe  strongly  that  modified  weak-shock  theory 
is  applicable  over  a  much  wider  range  of  experimental  conditions.   In 
particular,  the  theory  should  be  useful  for  source  sound  pressure  levels 
up  to  165  dB  and  should  not  have  a  lower  strength  limit. 

We  will  now  consider  our  noise  results.   Our  model  adequately 
describes  the  propagation  distortion  of  broadband  noise  waves  both  before 
and  after  the  formation  of  shocks.   The  frequency  band  limits  at  the  source 
in  our  experiments  were  700  Hz  and  2200  Hz  and  the  source  strength  limits 
were  155  dB  and  160  dB.   As  in  the  previous  case  we  believe  that  our 
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theory  is  applicable  over  a  -wider  range  of  experimental  conditions  than 
those  used  in  our  experiments. 

The  spectral  content  of  noise  waveforms  has  been  investigated 
theoretically  and  experimentally  for  various  propagation  distances.   Signif- 
icant spectral  shaping,  attributable  to  nonlinear  effects,  has  been  predic- 
ted and  these  theoretical  predictions  have  been  experimentally  confirmed. 
We  conclude, that  nonlinear  steepening  causes  a  shift  of  energy  to  higher 
frequencies  whereas  the  merging  of  shocks  and  the  attendant  waveform 
simplification  causes  a  shift  of  energy  to  lower  frequencies.   The  high 
frequency  build-up  commences  as  soon  as  the  nonlinear  wave  is  formed  at 
the  source.   The  low  frequency  build-up,  however,  is  associated  with  a 
more  mature  stage  of  the  distortion,  after  the  shocks  have  formed. 
Therefore,  we  expect  to  see  the  high  frequency  enhancement  at  considerably 
lower  SPL's  than  the  low  frequency  effect. 

Let  us  now  turn  to  the  body  of  material  included  here  that 
relates  to  Burgers  model  of  turbulence.  We  have  demonstrated  that 
use  of  weak-shock  theory  produces  several  results  that  are  in  agreement 
with  results  obtained  by  use  of  the  Burgers  turbulence  model.   However, t 
we  were  unable  to  verify  the  asymptotic  value  of  i   predicted  by  Tatsumi 
and  Kida  (1972).   The  general  agreement  between  the  sets  of  results 
indicates  that  the  two  theories  while  quite  different  in  outward  appearance 
and  originally  intended  application  are  actually  very  similar.   The 
failure  to  verify  the  predicted  value  of  |  merely  indicates  that  the 
value  £  =  l/2  cannot  be  verified  using  weak-shock  theory.   We  conclude 
that  if  I  approaches  l/2  it  does  so  at  wave  amplitudes  for  which  weak- 
shock  theory  is  invalid. 
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k.      The  Question  of  Dispersion 

In  Chapter  2  we  announced  our  intention  of  including  dispersion 
in  our  analysis.   Later  in  Chapter  3  we  indicated  in  a  more  complete  manner 
why  dispersion  was  important  in  our  particular  problem.   In  this  section 
we  shall  indicate  why  dispersion  is  important  in  some  instances  and  may 
be  ignored  in  others.   The  reader  is  cautioned  that  our  views  on  this 
subject  are  to  an  extent  speculative  and  are  not  as  yet  fully  verified  by 
laboratory  observations. 

Let  us  first  consider  finite-amplitude  waves  propagating  in  a 
pipe  that  are  too  weak  to  form  shocks.   A  wealth  of  experimental  results 
has  been  obtained  for  this  case.   The  experimental  work  includes  that  of 
Thuras,  Jenkins,  and  O'Neil  (1935),  Blackstock  (1962b),  Burns  (1966a), 
Pernet  and  Payne  (1969),  and  Cruikshank  (1966).   Each  of  the  above  research- 
ers either  ignored  dispersion  entirely  or  reported  that  its  effects  were 
slight.   In  work  where  source  levels  are  restricted  to  levels  insufficient 
to  form  shocks  only  the  first  few  harmonics  of  the  fundamental  wave 
frequency  reach  appreciable  levels.   Let  us  concentrate  our  attention  now 
on  the  second  harmonic.   The  observations  that  we  shall  make  about  the 
second  harmonic  can,  in  general,  be  extended  to  higher  harmonics.   One 
may  think  of  the  second  harmonic  as  being  generated  all  along  the  propaga- 
tion path  by  interaction  of  the  fundamental  with  itself.   This  concept 
is  often  referred  to  as  the  "virtual  source"  concept.   As  long  as  both 
fundamental  and  second  harmonic  have  the  same  propagation  speed,  the  second 
harmonic  signal  produced  at  a  distance  x,  say,  is  in  phase  with  all  the 
second  harmonic  signals  produced  at  distances  less  than  x  and  traveling 
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the  remaining  way  to  point  x  as  a  second  harmonic.   However,  if  there  is 
dispersion,  the  various  second  harmonic  signals  arriving  at  point  x  will 
not  be  perfectly  phased.  Now  in  the  case  of  dispersion  due  to  tube  wall 
effects,  the  dispersion  is  quite  small,  and  a  sizeable  distance  is 
required  for  the  imperfect  phasing  to  have  an  appreciable  effect.   In  the 
experiments  mentioned  earlier  the  distances  used  were  not  sufficient  for 
dispersion  to  be  very  noticeable.  For  weak  waves,  the  primary  signal 
components  are  the  fundamental  and  the  second  harmonic  (and  perhaps 
the  third  harmonic  as  well).   Since  the  dispersion  over  a  frequency  range 
of  only  a  few  octaves  is  very  small,  the  effects  of  dispersion  are  not 
noticeable  until  the  distance  becomes  great. 

Let  us  now  consider  waves  whose  source  levels  are  sufficient 
to  cause  the  formation  of  shocks.   In  this  case  the  fact  that  many  harmon- 
ics are  present  means  there  is  an  appreciable  spread  in  the  propagation 
speeds  of  the  wave  components.   Hence  we  expect  to  see  evidence  of 
dispersion  at  short  propagation  distances.   For  instance,  in  both  the 
work  reported  here  and  in  the  work  of  McKittrick  et  al.  (1967)  the  effects- 
of  dispersion  are  readily  apparent  in  the  form  of  waveform  asymmetry.   As 
source  levels  are  increased  beyond  approximately  155  dB,  however,  the 
signs  of  dispersion  close  to  the  source  disappear.   In  our  Fig.  (6-9), 
for  example,  the  waveforms  at  distances  less  than  6l  ft  do  not  exhibit  the 
rounded  peaks  characteristic  of  the  effects  of  dispersion.   The  reason 
is  that  the  nonlinear  steepening  effect  is  very  strong  and  it  keeps 
the  peaks  sharp  in  spite  of  dispersion.   In  other  words,  for  very 
strong  waves  dispersion  is  not  absent;  its  effects  are  merely  obscured 
by  the  more  powerful  nonlinear  effects. 
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Summarizing  the  above  discussion,  we  conclude  that  for  weak 
waves  and  relatively  short  distances  dispersion  may  be  ignored.   For 
waves  of  moderate  source  level  the  effects  of  dispersion  are  readily 
apparent  and  dispersion  must  be  included  in  the  analysis.   Finally,  for 
very  strong  waves,  the  effects  of  dispersion  are  small  compared  to  the 
nonlinear  effects  and,  at  short  distances,  can  probably  be  ignored. 

In  discussions  contained  in  this  section  we  have  been  purposely 
vague  about  the  sound  pressure  levels  and  distances  for  which  dispersion 
should  or  should  not  be  considered.   The  discussion  presented  here  is 
intended  as  a  general  set  of  guidelines  useful  in  the  analysis  of  propaga- 
tion of  finite-amplitude  waves  in  pipes.   A  more  definitive  study  leading 
to  a  classification  of  the  importance  of  dispersion  in  terms  of  source 
level,  frequency,  and  propagation  effects  needs  to  be  carried  out. 
5.   Avenues  of  Further  Experimental  and  Theoretical  Research 

As  is  usually  the  case  with  long  analyses,  we  came  to  the  end 
with  as  many  new  questions  as  we  have  found  answers  to  old  ones.   In  this 
section  we  shall  indicate  some  possible  avenues  of  continued  research 
that  may  build  upon  the  present  work. 

The  results  presented  here  are  limited  to  plane  waves.   In 
most  practical  propagation  problems,  however,  the  wave  motion  is  cylindri- 
cal or  spherical.   Therefore,  a  need  exists  to  modify  our  programs  so 
that  cylindrically  and  spherically  spreading  waves  may  be  considered. 
Since  weak-shock  theory  expressions  have  been  derived  for  both  cylindri- 
cally and  spherically  spreading  waves,  the  modification  should  not  be 
difficult  to  implement.   A  significant  simplification  will  accrue  in 
that  dispersion  is  negligible  for  waves  in  an  open  medium. 
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We  have  not  applied  our  theory  to  deterministic  waveforms  other 
than  sinusoids.   A  research  project  devoted  to  the  investigation  of  multi- 
ple frequency  source  waveforms  should  certainly  be  worthwhile.   The 
application  of  weak-shock  theory  to  waves  having  multiple  frequency 
components  at  the  source  is  very  difficult  to  handle  analytically. 
The  use  of  our  weak-shock  algorithm  will  greatly  reduce  the  time  required 
to  obtain  problem  solutions. 

In  Chapter  2  we  mentioned  Blackstock's  integro-differential 
equation.   Thus  far  there  are  no  known  solutions  of  this  equation.   In 
addition  it  is  not  known  whether  the  equation  could  be  used  alone  or 
whether  it  would  have  to  be  supplemented  when  shocks  form  by  use  of  the 
Rankine-Hugoniot  relations.   A  research  effort  aimed  at  obtaining  anlaytic 
or  numerical  solutions  to  Blackstock's  equation  might  be  fruitful. 

In  the  previous  section  we  discussed  in  very  general  terms  the 
importance  of  dispersion  in  various  pipe-borne  wave  propagation  problems. 
A  fuller,  more  detailed  experimental  and  theoretical  analysis  of  this 
topic  is  needed. 

Another  very  important  need  is  an  analysis  and  classification 
of  the  effects  reported  in  our  study  of  noise  waves.   This  analysis 
should  investigate  the  effect  of  varying  the  noise  bandwidth.  The 
information  obtained  should  be  reduced  so  that  engineering  estimates 
of  spectral  changes  can  be  made.   The  analyses  should  be  directed  also  at 
answering  whether  or  not  phase  information  is  important.   We  have  essen- 
tially thrown  phase  information  away.   Finally,  consider  a  case  where 
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we  have  a  noise  spectrum  but  not  the  noise  waveform.  How  good  a 
prediction  of  subsequent  spectra  may  be  made  based  upon  the  information 
given?  We  do  not  know  the  answer  to  this  question  as  yet. 

Lastly,  we  feel  that  the  work  begun  in  Appendix  E  needs 
attention.  We  have  indicated  that  for  very  strong  waves  ordinary  weak- 
shock  theory  is  inadequate.   However,  the  three-term  solution  and 
integral  solution  given  in  Appendix  E  cannot  as  yet  be  regarded  as  full 
theories  because  the  assumption  of  simple  flow  implicit  in  their  deriva- 
tion needs  to  be  investigated. 


APPENDIX  A 


LISTING  AND  FLOW  CHART  OF  SUBROUTINE  WAVEPROP 
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DRIVER  PROGRAM  LISTINGS  AND  FLOWCHARTS 

The  values  assigned  to  the  various  program  parameters 
in  the  cases  shown  are  example  values. 
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APPENDIX  D 


LISTING  AND  FLOWCHART  OF  PROGRAM  DOITALL3 


The  values  assigned  to  various  program  parameters 
in  this  case  are  example  values. 
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APPENDIX  E 

The  Limitation  on  the  Strength  of  the  Waves 
that  may  be  Considered  using  Weak-Shock  Theory 


Previous  researchers  [see,  for  example,  Rudnick  (1952,  Black- 
stock  (1964)  and  Temkin  (1969)]  have  used  a  variety  of  upper  bounds  for 
the  strength  of  the  waves  that  may  be  considered  using  weak-shock  theory. 
In  this  appendix  we  shall  analyze  the  limit  using  entropy  considerations, 
Our  analysis  will  basically  follow  Rudnick's  (1952).*   In  Rudnick's  deri- 
vation and  in  ours  the  implicit  assumption  of  simple  wave  flow  is  made. 
For  waves  of  high-intensity,  however,  reflections  from  the  shocks  may  be 
large  enough  to  invalidate  the  simple  wave  assumption,  even  as  an  approx- 
imation.  Therefore,  we  should  not  regard  the  results  obtained  here  as 
a  tool  for  calculating  the  attenuation  of  very  intense  waves,  but  rather 
as  an  estimate  useful  in  establishing  the  limits  of  weak-shock  theory. 
Our  results  may  be  used  in  problems  concerning  very  intense  waves  if 

1)  further  analysis  shows  that  the  reflected  waves  may  be  ignored,  or 

2)  they  are  used  in  conjunction  with  the  Riemann  invariants  [see,  for 
example,  Thompson  (1972)]. 

In  Rudnick's  derivation  of  the  weak-shock  propagation  of  saw- 
tooth waves  use  is  made  of  the  approximate  relation 


s2  "  si  =7T  6 


-  BR_  .3 

6Y  (E-l) 


*A  brief  but  more  accessible  form  of  this  derivation  appears  in  the  paper 
by  Rudnick  (1953).  2Qg 
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Equation  (E-l)  is  numbered  (3-17)  in  Chapter  III.   Use  of  Eq .  (E-l) 
leads  to  the  weak-shock  relation 


d(l/6)  _  y  +  1  (E-2) 

d(x/\)    2Y 


The  solution  of  Eq .  (E-2)  subject  to  the  boundary  condition  6=6   at 
x  =  0  is 


6  "  6     2Y    X  (E"3) 

o     ' 


If  we  use  the  exact  expression  for  the  entropy  change  across  the  shock, 
an  improvement  is  obtained.   The  exact  entropy  change  across  any  shock  in 
a  perfect  gas,  is  given  by 

As  =  Cv  in   [(1  +  6)(p2/Pl)"Y]     ,  (E"4) 

where  C   is  the  specific  heat  at  constant  volume.   The  absolute  density 
v 

ratio  p«/p,  is  obtained  from  Eq .  (3-13).   Use  of  the  perfect  gas  relations 
and  inversion  at  the  resulting  equation  yields  the  desired  relation 


P_2  =  6  (V  +  1)  +  2y  (E_5^ 

Pi   6(Y  -  1)  +  2y 


Equation  (E-4)  is  derived,  for  example,  by  Lindsay  (1960) 
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Following  Rudnick's   procedure  and  using  Eqs .  (E-4)  and  (E-5)  in  lieu 
of  Eq.  (E-l)  we  obtain  the  differential  equation 


d(l/5)  =  6_y_ 

d(x/x) 


in 


(Y  "  1)6' 


(1  +  6) 


6 (y  +  1)  +  2y 


6  (Y  -  D  +  2 


Yj 


(E-6) 


in  place  of  Eq .  (E-2) . 

We  have  obtained  approximate  solutions  of  Eq .  (E-6)  by  two 
different  methods.   In  the  first  method  the  right  hand  side  of  Eq.  (E-6) 
was  expanded  in  a  power  series  and  terms  through  the  fifth  order  in  6 
retained.   The  final  result  is 


d(l/6)  1  V  +  1   3(y  +  1)6    3(lly2  +  l)(y  +  1)62 
d(x/\)    2Y        4Y  40y3 


(E-7) 


Solution  of  Eq .  (E-7)  subject  to  the  boundary  condition  6=6   at 
x  =  0  yields  the  result 


T-T-'^ST1     f   -F(6)+F(6  ) 

6    6Q    2y     X  o 


(E-8) 


where 


F(6)  = 


r20Y 


30Y  (6  "  1) 


4y 


in 


+ 


+  3, 


132y  +  15 


2       1/2 
[1740y  +  240]i/z 


tan 


40  Y 


30  y 


2         1  /? 

,[1740y  +  240]  u- 
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Note  that  Eq .  (E-8)  is  an  implicit  result,  since  we  have  not  obtained 

an  explicit  solution  for  the  quantity  6(x). 

In  our  second  solution  Eq .  (E-6)  was  numerically  integrated  by 

means  of  Simpson's  rule.   In  Fig.  (E-l)  we  have  plotted  the  quantity 

1/6  versus  x/\  for  y   =  1.41  using  Eqs .  (E-3),  (E-8)  and  the  results  of  the 

numerical  integration.   Three  values  of  6  "  have  been  used,  6  "  =  2,  10 

o  o 

and  15.   For  each  value  of  6    the  top  curve  is  the.  ordinary  weak-shock 

result  Eq .  (E-2),  the  middle  curve  is  the  three-term  result  Eq .  (E-8"),  and 

the  bottom  curve  is  the  result  of  the  numerical  integration.   In  each  case 

the  three  curves  shown  appear  to  be  diverging.   In  reality  if  the  plots 

are  extended  to  larger  values  of  x/\  the  three  curves  become  parallel  with 

the  common  slope  given  by  the  quantity  (y  +  l)/2y- 

Note  that  the  bottom  curve  in  the  set  for  6=2  has  consider- 

o 

able  curvature  for  small  values  of  x/\.   This  curvature  indicates  that 
very  intense  waves  will  attenuate  at  a  rate  considerably  less  than  that 
predicted  by  use  of  weak-shock  theory. 

The  upper  limit  of  applicability  of  weak-shock  theory  was  ob- 
tained by  evaluating  the  right  sides  of  Eq .  (E-2)  and  (E-6)  as  functions  of 
6  and  calculating  the  percentage  difference  between  the  two  results.   The 
percentage  differences  along  with  results  obtained  by  evaluating  the  right 
sides  of  Eqs.  (E-2),  (E-6)  and  (E-7)  are  shown  in  Table  (E-l).   The  results 
obtained  using  Eq .  (E-7)  are  included  for  the  purpose  of  comparison  with 
the  results  obtained  using  Eq .  (E-6). 
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TABLE  E-l 

Analysis  of  the  Upper  Strength  Limit  of  Weak-Shock  Theory 

6  A  BCD 


Weak  Shock 
Theory 


Eq.  (E-7)    Eq.  (E-6)    7  Error 


0.5  0.855  0.582  0.7938  88.32 

0.2  0.855  0.658  0.647  32.14 

0.1  0.855  0.741  0.740  15.54 

0.0667  0.855  0.776  0.775  10.32 

0.05  0.855  0.7942  0.7938  7.71 

0.01  0.855  0.8419  0.8418  1.65 


Notes:   1)  y   =   1.41 

2)   D  =  (A  -  C)/C  X  100 

We  see  from  Table  E-l  that  the  error  involved  in  using  the  weak-shock 

result  when  6  =  0.01  is  less  than  2%   whereas  the  error  when  6  =0.1  is 

15.57o.   Clearly,  the  use  of  weak-shock  theory  will  produce  results  of 

doubtful  validity  for  values  of  6   larger  than  0.1.   The  value  6   =10 

o  o 

corresponds  to  a  source  SPL  of  about  -\jfo   dB.   For  waves  stronger  than  165  dB 
improved  versions  of  weak-shock  theory  based  upon  Eq .  (E-6)  in  conjunction 
with  the  Riemann  invariants  should  be  used. 
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